{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Random Survival Forests\n",
    "\n",
    "This notebook demonstrates how to use [Random Survival Forests](https://scikit-survival.readthedocs.io/en/latest/api/generated/sksurv.ensemble.RandomSurvivalForest.html#sksurv.ensemble.RandomSurvivalForest) introduced in [scikit-survival](https://github.com/sebp/scikit-survival) 0.11.\n",
    "\n",
    "As it's popular counterparts for classification and regression, a Random Survival Forest is an ensemble\n",
    "of tree-based learners. A Random Survival Forest ensures that individual trees are de-correlated by 1)\n",
    "building each tree on a different bootstrap sample of the original training data, and 2)\n",
    "at each node, only evaluate the split criterion for a randomly selected subset of\n",
    "features and thresholds. Predictions are formed by aggregating predictions of individual\n",
    "trees in the ensemble.\n",
    "\n",
    "To demonstrate Random Survival Forest, we are going to use data from the German Breast Cancer Study Group (GBSG-2) on the treatment of node-positive breast cancer patients. It contains data on 686 women\n",
    "and 8 prognostic factors:\n",
    "1. age,\n",
    "2. estrogen receptor (`estrec`),\n",
    "3. whether or not a hormonal therapy was administered (`horTh`),\n",
    "4. menopausal status (`menostat`),\n",
    "5. number of positive lymph nodes (`pnodes`),\n",
    "6. progesterone receptor (`progrec`),\n",
    "7. tumor size (`tsize`,\n",
    "8. tumor grade (`tgrade`).\n",
    "\n",
    "The goal is to predict recurrence-free survival time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from sklearn import set_config\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import OrdinalEncoder\n",
    "\n",
    "from sksurv.datasets import load_gbsg2\n",
    "from sksurv.ensemble import RandomSurvivalForest\n",
    "from sksurv.preprocessing import OneHotEncoder\n",
    "\n",
    "set_config(display=\"text\")  # displays text representation of estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to load the data and transform it into numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = load_gbsg2()\n",
    "\n",
    "grade_str = X.loc[:, \"tgrade\"].to_numpy(dtype=\"str\")[:, np.newaxis]\n",
    "grade_num = OrdinalEncoder(categories=[[\"I\", \"II\", \"III\"]]).fit_transform(grade_str)\n",
    "\n",
    "X_no_grade = X.drop(\"tgrade\", axis=1)\n",
    "Xt = OneHotEncoder().fit_transform(X_no_grade)\n",
    "Xt.loc[:, \"tgrade\"] = grade_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the data is split into 75% for training and 25% for testing, so we can determine\n",
    "how well our model generalizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_state = 20\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(Xt, y, test_size=0.25, random_state=random_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training\n",
    "\n",
    "Several split criterion have been proposed in the past, but the most widespread one is based\n",
    "on the log-rank test, which you probably know from comparing survival curves among two or more\n",
    "groups. Using the training data, we fit a Random Survival Forest comprising 1000 trees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomSurvivalForest(min_samples_leaf=15, min_samples_split=10,\n",
       "                     n_estimators=1000, n_jobs=-1, random_state=20)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rsf = RandomSurvivalForest(\n",
    "    n_estimators=1000, min_samples_split=10, min_samples_leaf=15, n_jobs=-1, random_state=random_state\n",
    ")\n",
    "rsf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check how well the model performs by evaluating it on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.67453'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_index = rsf.score(X_test, y_test)\n",
    "f\"{c_index:.5f}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives a concordance index of 0.68, which is a good a value and matches the results\n",
    "reported in the [Random Survival Forests paper](https://projecteuclid.org/euclid.aoas/1223908043)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting\n",
    "\n",
    "For prediction, a sample is dropped down each tree in the forest until it reaches a terminal node.\n",
    "Data in each terminal is used to non-parametrically estimate the survival and cumulative hazard\n",
    "function using the Kaplan-Meier and Nelson-Aalen estimator, respectively. In addition, a risk score\n",
    "can be computed that represents the expected number of events for one particular terminal node.\n",
    "The ensemble prediction is simply the average across all trees in the forest.\n",
    "\n",
    "Let's first select a couple of patients from the test data\n",
    "according to the number of positive lymph nodes and age."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>estrec</th>\n",
       "      <th>horTh=yes</th>\n",
       "      <th>menostat=Post</th>\n",
       "      <th>pnodes</th>\n",
       "      <th>progrec</th>\n",
       "      <th>tsize</th>\n",
       "      <th>tgrade</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>119</th>\n",
       "      <td>33.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>574</th>\n",
       "      <td>34.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>421</th>\n",
       "      <td>36.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>76.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>65.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>80.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226</th>\n",
       "      <td>72.0</td>\n",
       "      <td>1091.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      age  estrec  horTh=yes  menostat=Post  pnodes  progrec  tsize  tgrade\n",
       "119  33.0     0.0        0.0            0.0     1.0     26.0   35.0     2.0\n",
       "574  34.0    37.0        0.0            0.0     1.0      0.0   40.0     2.0\n",
       "421  36.0    14.0        0.0            0.0     1.0     76.0   36.0     1.0\n",
       "24   65.0    64.0        0.0            1.0    26.0      2.0   70.0     2.0\n",
       "8    80.0    59.0        0.0            1.0    30.0      0.0   39.0     1.0\n",
       "226  72.0  1091.0        1.0            1.0    36.0      2.0   34.0     2.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_test_sorted = X_test.sort_values(by=[\"pnodes\", \"age\"])\n",
    "X_test_sel = pd.concat((X_test_sorted.head(3), X_test_sorted.tail(3)))\n",
    "\n",
    "X_test_sel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predicted risk scores indicate that risk for the last three patients is quite\n",
    "a bit higher than that of the first three patients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     92.485064\n",
       "1    101.425923\n",
       "2     73.756955\n",
       "3    173.702023\n",
       "4    172.987572\n",
       "5    154.022612\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(rsf.predict(X_test_sel))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a more detailed insight by considering the predicted survival function. It shows that the biggest difference occurs roughly within the first 750 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAELCAYAAADURYGZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAABDTElEQVR4nO3de1zUZdr48c8IKh7AE4ICnghRQNAS1DQNLdOIdNMUs7R0jeywafXs5m8Pmdaz1W7tVmu7Rge3Nk2zzYdSozKl8pCIZnkEVDRBBBVHQESYYX5/jN9pZphhBpgjc71fr14539PcN6NzcZ+uW6XT6XQIIYTwWW3cXQAhhBDuJYFACCF8nAQCIYTwcRIIhBDCx0kgEEIIHyeBQAghfJy/uwvQVMHBwfTv37/J912+fJlOnTo5vkAexhfq6Qt1BN+opy/UETyjnidPnuT8+fMWz3ldIOjfvz+5ublNvi87O5vk5GTHF8jD+EI9faGO4Bv19IU6gmfUMzEx0eo56RoSQggfJ4FACCF8nAQCIYTwcV43RiCEEO5SV1dHUVERNTU1TbqvS5cuHDlyxEmlMhUQEEBERARt27a1+x6nBYL58+ezceNGQkJCOHjwYIPzOp2ORYsWsXnzZjp27Mi///1vbrjhBmcVRwghWqyoqIjAwED69++PSqWy+77KykoCAwOdWDI9nU7HhQsXKCoqYsCAAXbf57SuoQceeICsrCyr5z///HMKCgooKCggIyODhx9+2FlFEUIIh6ipqaFHjx5NCgKupFKp6NGjR5NbLE4LBOPGjaN79+5Wz2dmZjJ37lxUKhWjRo1CrVZTUlLirOIIIYRDeGoQUDSnfG4bIyguLqZPnz6G1xERERQXF9O7d2+nvN+Bt/7D/pWrGr1G21HD7lGHAKho05WLfj0I7tSekKD2Da5NiUxhRvQMp5RVCCGsycrKYtGiRWi1WhYsWMCSJUta/Ey3BQJL++FYi2QZGRlkZGQA+j667Ozs5rxho6e1ugtwGUZvHQro0HWoY/vIw3AVqi6aXnu6XT3VZSfoeaZn08vhZFVVVc37+XgRX6gj+EY9va2OXbp0obKyssn3abXaZt1n6TkPP/wwmZmZhIeHk5yczC233MLgwYNNrqupqWnSz9VtgSAiIoLTp08bXhcVFREWFmbx2vT0dNLT0wH96rjmrtBr7L5Pl7zIqZ/zAajVqml3JZB/1gZzvuqqyXWVNRr+GaZG27bS7SsFLfGEFYzO5gt1BN+op7fV8ciRI80a9HXUYPGuXbuIjo4mISEBgNmzZ7NlyxaSkpJMrgsICOD666+3+7luCwRTpkxhxYoVzJo1i927d9OlSxendQvZVZ4Xf2le/WP2fOq0F8m+mETMmGQSbp1sOLdm989c/uFOSlRXmbf+dlKGzpcuIiGES1jqUt+9e3eLn+u0QHDPPfeQnZ3N+fPniYiIYNmyZdTV1QGwcOFCUlJS2Lx5M1FRUXTs2JFVqxrvv3elDm37o8GfosMHKTp8kCM7sg3n/IBBdYMIiviB3OoicnctZ/POl0gJG8uM2/7utjILIVxr2WeHOHymwq5rtVotfn5+Nq+LDQti6Z1xVs83pUu9KZwWCD788MNGz6tUKt544w1nvX2LRLdpw2nVCM53rEJbe5Tzp6sM52qrzxLapTf9e8wi9vLn7AnUkNe2Bs58h7QLhBDO1JQu9aaQlcUWxE+Opt/GjRw7G0RZxDgCjAZiTu1/C015IaPiH+N/j6dw5EwFMWFPU+tX78YSCyFcrbHf3M05aowgKSmJgoICCgsLCQ8PZ+3ataxZs6bFz5VAYEG3tJl0S5sJc+YSVZNFv6dmG8699/T1nD95mq/eWsHEW2dD7wFo63X46aph1R36i+LvhsR5biq9EKK18vf3Z8WKFUyaNAmtVsv8+fOJi7M/IFl9rgPK1qrVHD3KqTlzDa8H+Eej7ngrmuotBJ75iXVLZzP57W6c9rvAPFUp1F6Gg2+QEtRZBpGFEA6XkpJCSkqKQ58p2UcbEZSaatItBBB68FNC/ENQ+UdwWV0LQD23U301ksO6flymPbmqqyzftZz1+evdUWwhhGgSaRE0wtBFZOTUnLmEa05QAlyp1AeCB+LvIXP/OKiF2JJPmBG4muXBPdh8YrO0CoQQHk9aBM3Qc9dq2vnrqKk8Sda/PmT2yL6se+hG1j10I4d7TyP2aiSJV2qg8qy7iyqEEDZJIGiioNRUAHqW6ZeLH8pezU9bTLOsfnB5BAB5l4uZlzWPeVnzpJtICOGxJBA0Ube0mfRatoyEwh/oEqD/ws/d9KXh/NRh4XyovYVhNd0YRDv9+dJcGTMQQngsCQTN0C1tJh2Tkohu0waVf4RhrABg9si+jBzQnVur27GqpIxVJWU8U69Px735x3fdVWQhhLBKAkELGMYKqkp4Y8Fjhv+uy8vlg8sjKO00EIAZdCbxqobc6iJpFQghWmT+/PmEhIQwZMgQhz1TAkEzGY8VtG8fbDheU1VCcOkhDlXfxuPtn4d5m2DeJlLahgCw+cRmt5RXCNE62Nr9sTlk+mgzGaaVLl1Kx+Ak+r39JgDvPvEk6rJqRvoFsLakgrQ3dwHwzAUtCaFtySvPY17WL6uOZYMbIURTjBs3jpMnTzr0mdIiaAFlrMBYp67t0GmK6FJ3iNjeQYbj1bVafqU+y6C2XQzHZBBZCOEJpEXgAMZpKEK0NRQBtecP8Ezsr4gbGw5Axt9vJf3S64wrqSV0hj7l9vr89SzftVwWngnhjT5fAmcP2HVpB60G/Oz4uu0VD7e/2MKCNZ0EghZSxgoUYQUnKYi6jkrtOb5Y+RzZ7+l/xG36X8/39THEXzyiT04XfzczEufJmIEQwu0kELSQeRqKU3Pm0u/4UU7HDuXqtR9vTVUJASfhp7iJ+F3aQtKp7XBqOxz4GDrXQ2AvN5VeCNFsTfjN/YqD0lA7iwQCBwtKTaXv0j0MvnqVfv95C4A3FjwGQOfRD/Ly/hRiSz7h7na7iDu1HXqFkFdzTgaQhRB2sbT7469//esWPVMCgYN1S5tJxcaNJuMG9Vp/arVldN+3g3UP30PamzCr5FZ+02U7KZffgYBfBpDzyvMAJBAIISyytftjc8isIScwT18dUqUF4PjeHYA+DUVs7yD+XDqKPpf68mJFF1ZNXsWqyasY1H0QuaW5MpNICOEy0iJwggbpq+fM5VNtiOHl7JF9mT2yL2t2/wyb4HzVVUKvnUuJTDFMKzUeSJbuIiGEs0ggcKPZI/ty6GvTj0D5sjcOArmlueSW5hqOSVAQQjiSBAIX0Wk16OrqLJ7rX3fCMKWUxHnMiJ5h8kW/Pn+9IQgoQQFkHEEI4RgyRuACyloDS4FgR4fx+i0uf/6B0p0fWLx/RvQMwxjCMzc+A8DyXctlnwMhhENIIHCBbmkzUVlZVdh59IO83PtvHNb1I7Q8F3JXNfqsGdEzeObGZ0gMTSSvPE8WpAkhWkwCgQvV1l/i3SeeZN2yJYZdzZRtLn/ocqv+oo2L7QoGMsNICN9z+vRpxo8fT0xMDHFxcbz22msOea4EAhfp7dcJlV9P1GXVnMk7ZrKrGehbBv+vTr8oxFoXkbmUyBRAUlsL4Sv8/f155ZVXOHLkCN9//z1vvPEGhw8fbvFzJRC4SBKXSGwXwYAb0tGpglGXVbPhlX1seGUfh74rZvbIvsRPWcz39TGcr7pq1zNnRM8gMTRRWgVC+IjevXtzww03ABAYGEhMTAzFxcUtfq4EAhfqp8nnrqduoGtIR1S685w58i6F+zL4YuVzZP3rQ2aP7EtggD9xtQdsdg8plFaBpLMWwrecPHmSH374gZEjR7b4WTJ91A0S77iNIzuyAbisruXimXwOZa8mbGA3vu8wXh8IDnwMifMafxC/TCFVFqDdH3C/M4suhLjmpZyXOFp+1K5rtVotfn5+Nq8b3H0wT4942uZ1VVVVTJ8+nVdffZWgoCCb19sigcCFqvfs4eK6j0hIm0nCrZMNx997OoPzJz/lyI5svu41lTFXthHXhOfOiJ4h4wRC+Ii6ujqmT5/Ovffey7Rp0xzyTAkELhKUmkr1nj2cXbqUio0bCUpNNaShCApJ5ELRPooOH6Sntg8EoN/wYtUdvzzg2mKzxuSV5/Fam9d4L+s9QFYgC+FM9vzmrqh0UBpqnU7Hr3/9a2JiYnjyySdb/DyFBAIXUb70KzZupHrPHqr37KFi40YAIm6Yzs/tBqPRFNG76CAf970R2gEllwju3F6/vuDUdv2DrAQDZaxArVYDsgJZiNZox44d/Oc//yE+Pp5hw4YB8Oc//5mUlJQWPVcCgQspyegurvvIEARqjh4lhP9ya/oLfLHyKB11Kr7vPY3lTGN3YTlUwO9Dvyf90uuU7vyAUCuBQElLkZ2dTXJysmyDKUQrdNNNN6HT6Rz+XAkEbmCcnfTUnLlU79lDWOoO2nfQfxzrHroRgDW7fyZzfzFfk0LCxa/g/GUef3MXU4eFM3tk30bfQ8YNhBD2kkDgZsZjB7qR41C1bWs4p6SrBih9vROdLx7hieInoBh2/zyDkTOeclexhRCtiKwjcLNuaTPptWwZYDkpnSJ09H106ns9A4I7Eas6ReeCDa4qohCilZNA4AG6pc2kY1ISOq2GmsqThjxEJhLnwbxNhD7+NafbX2f3s2XVsRDCFqcGgqysLAYNGkRUVBQvvvhig/OXLl3izjvvZOjQocTFxbFqlX2raVujoNRUOhAOwFdvrbAcDJpIVh0LIezhtECg1Wp59NFH+fzzzzl8+DAffvhhg+RIb7zxBrGxsfz4449kZ2fz1FNPUVtb66wiebRuaTPpHBRPl4ARgD4YGGcpNVddqyXtzV367S6tUFJWgySmE0JY57RAkJOTQ1RUFJGRkbRr145Zs2aRmZlpco1KpaKyshKdTkdVVRXdu3fH39+3x687tYti4oOPERE7hHMnCw2pKIwFd25Px3Z+7C4s5/cbDtgMBsreBbKRjRDeraamhhEjRhh6UZYuXeqQ5zotEBQXF9OnTx/D64iIiAZZ8h577DGOHDlCWFgY8fHxvPbaa7Rp49vDFuV+vfBrH0/a0hfp2X+AxWtCAwOIqz3AuuH6PCeZ+xvPPpgSmWLYu0C6iYTwXu3bt2fr1q38+OOP7N+/n6ysLL7//vsWP9dpv35bWvSgUqlMXn/xxRcMGzaMrVu3cvz4cSZOnMjYsWMbJFHKyMggIyMDgKKiIrKzs5tcnqqqqmbd50rdKw5S3q0XuetzOKftZlglbF7u3u0TGMR2BhV9wqBuf0StVhuusVTPnvTk/oD7ieoexdrytSzftZw1+9aQ2CmRMYFjnF8xB/OGz9IRfKGe3lbHLl26UFlZ2eT7tFpts+6zprKykurqaq5evUp1dXWDZ9fU1DTp5+q0QBAREcHp06cNr4uKiggLCzO5ZtWqVSxZsgSVSkVUVBQDBgzg6NGjjBgxwuS69PR00tPTAUhMTCQ5ObnJ5VFW3Hqyi6VlnN+Qj7prND39BlLatSuAhXInw6qf6HpqO9O77OTrjikkJ+sXoTVWz2SSic6PZvOJzeSW5nLs6jGOtT3mdTmJvOGzdARfqKe31fHIkSPNyhnkqFxDoA8qw4cP59ixYzz66KNMmDChwTUBAQFcf/31dj/TaYEgKSmJgoICCgsLCQ8PZ+3ataxZs8bkmr59+/L1118zduxYSktLycvLIzIy0llF8njd0mbSJ+t51ESTvTqPDh1q6dS1neWL4++GU9tJv/T6tQM32vUeSiqK9fnr2XxiM3nleYbjQgj7nf3zn7l6xL401BqtlnI70lC3jxlMr9//vtFr/Pz82L9/P2q1mrvuuouDBw8yZMgQu8phjdM65P39/VmxYgWTJk0iJiaGmTNnEhcXx8qVK1m5ciUAf/rTn9i5cyfx8fHccsstvPTSSwQHBzurSF4hfnI0g/L0AVNdVk3R4YPW1xWkvgrA9Ze22JxBZM5432MhhPfp2rUrycnJZGW1fKq5U6fopKSkNMiKt3DhQsOfw8LC+PLLL81v82nd0mYyHGiXtZMfrmUkPbIj22T/AoPEeZTu/ICOVVc5XFIBwMPyvS6ES9j6zd2Yo7qGzp07R9u2benatStXrlxhy5YtPP20/emwrfHtKToeqlvaTMat+iMh/iGo/COstwq4NoNI9TNr2z3HLdWyVkCI1qykpITx48eTkJBAUlISEydOJDU1tcXP9e1J+x4uXHOCMlutgvi7Aej/8w9U12r54+7R/CvPvgylCmWNgcLbBo+F8BUJCQn88MMPDn+uzRbBwYMHHf6mwj79NPkEXwlA5R/BZbWVFdfXchBVdYuhYzv9YJSy0MyecQNljYFCWWsgi8+E8B02WwQLFy6ktraWBx54gNmzZ9P12pRG4XxBqamErvyS0jB/rlQ2nnojNDCA0PJcXo/+lr3JD5O5v9gwbtBYy0CZRaSQ2URC+B6bLYLt27ezevVqTp8+TWJiIrNnz+arr75yRdl8Xre0mQyMqMVfU91oimrA0EUUWvots0f2Zd1DNxLbO6jxeyyQ2URC+B67BosHDhzI888/z0svvcQ333zD448/zuDBg/nkk0+cXT6fF3RtIOhqTXHjGUkT50G/mxz63pLCWgjfYLNr6KeffmLVqlVs2rSJiRMn8tlnn3HDDTdw5swZbrzxRqZNm+aKcvqsbmkz6bj5KJdqcqwPGDtBSmSKYbzAUuZSGVAWovWwGQgee+wxHnzwQf785z/ToUMHw/GwsDCef/55pxZO6HVqF0WF5oz1AWMjnasKYdUdANxSfQNfd0yxcYdlype8pSCQW5pLbmkum09sloAgRCtgMxBMmzaNOXPmmBx77bXXWLRoUYPjwjnCNScoAdRnT7Ju2RJixiRbnUpapVbTFeDsAcboLvHn0lGs2f2z3VNJjZkPJCuUAWUlICjXCiFcQ6vVkpiYSHh4OBs3bmzx82yOEbz//vsNjv373/9u8RsL+8VPjiZQ0422qq4UHT5ofQezxHnsv/5/Yd4m6BVPXO0B7vH72maa6qZSBpRl0xsh3OO1114jJibGYc+zGgg+/PBD7rzzTgoLC5kyZYrhv/Hjx9OjRw+HFUDYpuxeFtzpViY++BiAxQ1rTFybRXRfpxx2F5Y3KQ+RvWTTGyFcr6ioiE2bNrFgwQKHPdNq19Do0aPp3bs358+f56mnnjIcDwwMJCEhwWEFEE2TcOtkjuzINqSdsDp4nDgPDnxMcGUNVOg3r2lO95Atyr7IxuMGynHpLhLC8RYvXsxf/vIXh+5vYDUQ9OvXj379+rFr1y6HvZlomYo23dnwyj66hAyl6PBBu2YRhQYGMHJAd0OrwNHBwDytNUhQEL7hu4/yOX+6yq5rtVotfnakoQ7u05mxM6Otnt+4cSMhISEMHz7coRv6WA0EN910E9u3bycwMNBkZzGdTodKpaKiosJhhRC2hWtOUF9bzZmCdkAoXQJ7233v1GHhhrQT0PhK4+YyHlg2DgqyQlkIx9mxYweffvopmzdvpqamhoqKCu677z4++OCDFj3XaiDYvn07gEObH6L54idH02/jRk5dPcuB9qO5Uq3F3qS2yhf/7zcccFoXkTHjoGCczE6I1qSx39zNOSoN9QsvvMALL7wA6Hd3e/nll1scBKCRQFBeXt7ojd27d2/xmwv7dUubSbe0mfQDih9cw3mwPU5gZPbIvob8Q2lv6rv7mpKhtCUku6kQns1qIBg+fDgqlcrqJvQnTpxwasFE4+rbJ0D1FtvjBKe2Q+4qSJzH1GHhhsP2JKRzBGUwWSFdRUI4RnJyssP2e7YaCAoLCx3yBsLxwjUnKG8/Gn+/Y5w7Wci6ZUsAiBmTDP4Bv1x4bV9jDnwMifOYPbKv4Ys/7c1dhtaBM1sG5ovSpKtICM9jNRAcPXqUwYMHs2/fPovnb7jhBqcVSjSunyaf01eCqe0+gI6hnQF9N1HR4YP0vXkiKL8lXJtCatwqUCitA1e1DIQQnstqIPjb3/5GRkaGyRoChUqlYuvWrU4tmLAuKDUVNpQTUNebtKX6z+enLVl89dYKyguOmF5s1ipQKK2DtDd3OW1qqRDCO1gNBBkZGQBs27bNZYUR9umWNpM2W9ZQ7teLQ98VEzc23LDQTK1Wm16stArOHtAno4u/u0HLQJlaqqSicNUgshDCM9hMOldTU8M///lPtm/fjkqlYuzYsSxcuJCAgABbtwonCtecoNyvF/k5pcSNDW/84mvpJjirX0dg3jIADEFgd2E5uwvLydxfLAFBCB9hMxDMnTuXwMBAfvOb3wD6HERz5sxh/XrJK+NO/TT5FPtHcr7Inw2v6MdxGt3XOHGevkVgYbzAeBB5ze6fydxfbAgIynlHUja8kZlDQngGm4EgLy+PH3/80fB6/PjxDB061KmFEvbpeWwr/iNnAJ05X1RFbWUtdVfPmcwiMplaqowXbFys7y4y6yaCX4LCmt0/8/sNB0y6jIw1t7WgbHiz+cRmCQRCNEP//v0JDAzEz88Pf39/cnNzW/xMm4Hg+uuv5/vvv2fUqFEA7N69mzFjxrT4jUXLBKWmEr50KeH/t4OOSUnsDJhMmTaKtp00wC+ziIBfgoHypa/MJFIGka0EBMBiEGhJ99GM6BlsPrHZZJGZLDATomm2bdtGcHCww55nNRDEx8ejUqmoq6vj/fffp2/fvqhUKk6dOkVsbKzDCiCap1vaTAAqrm1KEVK0k/J+0+jYM4G0pRMMs4gaLDhTuolyV/0yiKwcN2PcZWSspd1HxovMzBPUWbpWgoQQzmU1EDhi1xvhXEraCQDmzKVMe5YaegG/pKu2ynjcQJlRpLDQQjBmrfvI3taBtQR15mQVshANqVQqbrvtNlQqFQ899BDp6ektfmajaaiNlZWVUVNT0+I3FM5TX1lBG6N9pcGOfETKjCKFjS4jY8bdR81dmGZtO0zQr0KWgWXhqbb9O4OyU/al2tFqtPj5205DHdIvkvEPNP7FvmPHDsLCwigrK2PixIkMHjyYcePG2VUOa2yOEXz66ac89dRTnDlzhpCQEE6dOkVMTAyHDh1q0RsLx1IWmbUxyhYbMybZ9r4FSstAoXQZKQFBucYK44VpxgntjLV0YHn5ruVsPrFZuomEAMLCwgAICQnhrrvuIicnx/mB4E9/+hPff/89t956Kz/88APbtm3jww8/bNGbCsdTFplVqrqx4ZV9RI8Itd09ZInxGMLGxb/MMDJmoaVgnNDOWEtSWChf+ptPbDYZS1Cr1ZzLPydBQbiVrd/cjTkqDfXly5epr68nMDCQy5cv8+WXX/LMM8+0+Lk2A0Hbtm3p0aMH9fX11NfXM378eJ5++ukWv7FwvHDNCTT05UxBe84UqA3Hm5Ku2sB4hpExKy0FawPLLU1uZ2kHtGNXjxlaCbZIK0K0JqWlpdx1110AaDQaZs+ezeTJTfh3bYXNQNC1a1eqqqoYO3Ys9957LyEhIfj727xNuEE/TT5d1Dn4/+5fZK/OIz+n1L7uIWvMu43gl5aCWe4ia5SWQktXLBuPJfzvZ//LsbbHbN6jtCKU+4XwdpGRkSbruhzF5jd6ZmYmAQEBvPrqq6xevZpLly45pCkinMO/qIjOGf+PkPB7OVOgJnpEPBGxQwzpqhssMmsqG7mLzBnPMGrJoLKxMYFj+EPyH2xetz5/vaHlIIFACOtsBoJOnTpx9uxZcnJy6N69O5MmTaJHjx6uKJtooqDUVNRqNdV79tCtqB1lg2YbWgUA507q95hoUSCAX2Ya2TmgDKaDyq6iLF4TQjTOZiB4++23Wb58ORMmTECn0/Gb3/yGZ555hvnz57uifKIJuqXN5MfQEIaWlsHSpZSGJqI5N4CEWyeTcOtk1i1b0rzxAnPNHFBWuDrttflWmQoZPxBCz2Yg+Otf/8oPP/xgaAVcuHCB0aNHSyDwYIZFZhvK0Vy4YDjeovECS2wNKCvHjYKCq9Nem2+VqZDxA9FcOp0OlUrl7mJYZWl7YVtsBoKIiAiTaU+BgYH06dPHrodnZWWxaNEitFotCxYsYMmSJQ2uyc7OZvHixdTV1REcHMw333zThOILaxrbs8AhrQKFtQFlJQiYpbBoLO21JS0NEtYWrMn4gWiOgIAALly4QI8ePTwyGOh0Oi5cuNDkbQIa3aEMIDw8nJEjRzJ16lRUKhWZmZmMGDHC5oO1Wi2PPvooX331FRERESQlJTFlyhSTPEVqtZpHHnmErKws+vbtS1lZWZMKLxpnac8CpVXw1VsrAAeMF1hiHBwspLCYHX83sx/Sn1cGkS1x5jaayviBrFwWTREREUFRURHnzp1r0n01NTUu28MlICCAiIiIJt1jNRBUXluhet1113HdddcZjk+dOtWuB+fk5BAVFUVkZCQAs2bNIjMz0yQQrFmzhmnTptG3r/4fekhISJMKLxqn7FkAXQ3HlC9+iwnpnME8hYWFFoK1L3pnDyxLSmzRVG3btmXAgAFNvi87O5vrr7/eCSVyDKuBYOnSpSavKysrUalUdO7c2a4HFxcXm3QhRUREsHv3bpNr8vPzqaurIzk5mcrKShYtWsTcuXObUn5hQ31lBRqz316ateK4ucy7joyT27mZeUpsGTwWvsrmGMHBgweZM2cO5eX6dMPBwcG8//77xMXFNXqfpQEL8z41jUbD3r17+frrr7ly5Qo33ngjo0aNIjo62uS6jIwMwx7KRUVFZGdn2yp2A1VVVc26z9sY17ND9EA4BDWlpQ3qrlaruXL+HG8+sRCA7gNj6Bnr/A2HhqnVdL10kLw1/4+SsEmNXqtWXwFoUHZHfpZRdVGo26gNg8dr9q2xeF1ip0TGBLp2Hw5f+DvrC3UEz6+nzUCQnp7O3/72N8aPHw/o/1E++OCD7Ny5s9H7IiIiOH36tOF1UVGRIVmS8TXBwcF06tSJTp06MW7cOH788ccGgSA9Pd2QajUxMZHk5GS7KmcsOzu7Wfd5G5N6Jidz4sE1qP160dNvoMnext01NYZWwbmThdSfKyE5eZHzC9h5AWxczKD8fzLo6k/6Y1ammv4rbxe7C8s50yHSpPvIkZ9lMvrn2EqHfaztMbsWsTmSL/yd9YU6gufX02YguHz5siEIACQnJ3P58mWbD05KSqKgoIDCwkLCw8NZu3Yta9aY/rY1depUHnvsMTQaDbW1tezevZsnnniiGdUQ1igDxkrKiegRoYYZRMr4gLK1pUuYTzk1n2pq5MnuE0grHEzm/mKnrzmwlQ5biNbMZiCIjIzkueeeY86cOQB88MEHdg2W+Pv7s2LFCiZNmoRWq2X+/PnExcWxcuVKABYuXEhMTAyTJ08mISGBNm3asGDBAoYMGdLCKglj8ZOjqV25htLQRMry+6M5d86kZaBw6JRSW4zHDYynmho7tZ2RbOf3oY/zj5KbTAaO1eornOngugVpgGytKVo1m4Hg3XffZenSpUybNg2AcePGsWrVKrsenpKSQkqK6YKehQsXmrz+7W9/y29/+1t7yyuaqFvaTIaj39Jye5U/ZapIw7oChcMXmjWFpXUIYFi1PNVvJ1/3Nv07lHex3mRBWmMcsVjN2taaEhBEa9FoINBqtcyYMYMtW7a4qjzCCZQtLU/Nex41kSbrCsDFs4jsdS25XSiwbt6NJqee/c9XHKm2PXvNUesQLG2tKdtoitak0UDg5+dHx44duXTpEl26dHFVmYST9NPkc/pKMJpzTZ8H7UmS+7Tl2eQbbV6X9uYuh+c1UoKCjBuI1sRm11BAQADx8fFMnDiRTp06GY6//vrrTi2YcDxlO8syddsG3UOA41JVO9Kp7fpuIjv2PjBnKa+RctwRgcFSMjvpLhLeyGYguOOOO7jjDs9ZBCSar1vaTPpkPY+aaLJX67s2jFNPgANTVTtC/N2/zChqRiAwz2sEDXMbtWQ/ZXPSXSS8lc1AcP/991NbW8vRo0dRqVQMGjSIdu3auaJswgn6afKpO3uWo/2mcXjTQUMgME5VrbQMAPe2Dsw3wblmmFoNhV0t32O2JsE8hYVxbqOW7qds/oU/L2teg1aCtBCEN7AZCDZv3sxDDz3Eddddh06no7CwkDfffJPbb7/dFeUTDhaUmkrUxo2crTxBuS6CDa/sAzCsL1BaBqCfUqrMJgI3BQXzXEWNsbQmoZHA4OgxBPNWgrQQhLewGQiefPJJtm3bRlRUFADHjx/njjvukEDgpYxnEBXXdwRCOF9UBdBgodlPW7IMQUAJCuDibiML00v3W1ulab4mwSzBnTllDMFRC9bMWwmWWggKaSkIT2IzEISEhBiCAOgXmEmWUO/XT5Ov/++p2YZWgTnzoPDVWyucm766pZqY4G72yL6GfZTT3tzl8A1yrG2KIy0F4WlsBoK4uDhSUlKYOXMmKpWK9evXk5SUxCeffAJgWGgmvE/N0aOcmjOXmoDJJhvYWOLy9NWOYmPW0dRh+vo6Y+8Da2krzFsK0joQ7mYzENTU1BAaGmrYOaxnz56Ul5fz2WefoVKpJBB4qaDUVMOfQ4p2Ut5vWoOFZuacssOZMymzjsz3VDYaN1DGDJy994ExSyuVQVoIwn1sBgJ700kI76KMFQAwZy5l2rMYb2BjjVvTUTSVpT2VjQeUrWQ9dTbzlcqyZaZwN5uBQPiO80VVbHhln2EGkSUemY6iMebjBsqAshIQlGvAMFZgiaPHDxTK5jhCuFMbdxdAeIaex7bStf0VzhdVkZ9T6u7iOE/iPJi3CVJf1b++1lqYOiyc2N5BFm85XFJhV4I7IbyVtAgEQamphC9dysDiWvYNW2xoGQCNtg68mtlitdnA7OGWu4rS3txlsbXgyFaCpWmmarWa97Les3mvDDaLlrIaCP72t781euOTTz7p8MII9+iWNpOKjRsB/Re/wnh9gTmPzEvUVMaL1RpZc6DMLDLmyFlG1qaZ2kOmogpHsBoIKisrXVkO4SHixoYbvvg3vLLPYuvAI/MSNYfx+MGqO35JZWEjTQXg0FlG1qaZ2rO9YWOL1sxJy0FYYzUQLF261JXlEB5AWVcQlJpKt7SZJq2DMwVqzhSoAdO8RK2G0jqwsRrZmKNTXDeHva0JmaYqGmPXOoJ33nmHQ4cOUVNTYzj+7rvvOrVgwrWUdQU1R48C+u4i49bBoe+KyV6d1yBraauhtA6MWwbGzFoJllJcO2tmUWMa22vZmDJNVZmqKq0DYcxmIJgzZw6DBw/miy++4JlnnmH16tXExMS4omzChQw5iObMtXhe+eI3Dgbg4r2OXcFSkjsLrQTzFNdKemvjc55E+dLffGKzyXabziTBxnvYDATHjh1j/fr1ZGZmcv/99zN79mwmTZrkirIJN1G6iMyFpaaSfO8YslfnkZ9TalhcpqSd8OqBY4WlPZSttBJmx9/N7If0167Z/TO/33DAaZvgOILSelC223QmGcT2LjYDQdu2bQHo2rUrBw8epFevXpw8edLZ5RJuYpx6wlj1nj1U79lD56QkQsLvBX4ZJFbSTigrjltFQDBmqZVgtiDN1iY4nhgQnEm28vQuNgNBeno6Fy9e5LnnnmPKlClUVVXx3HPPuaJswg1MUk8YubjuIyo2bqR6zx5q60YbjisDx0rKarelq3YmS62E3FX6HEZGeYzM1yIom+A4I6GdN5AU3N7DZiCYN28efn5+3HzzzZw4ccIVZRIeSAkQF9d9BBvKKb+oM0lHYRwQvC5DaXPYymPEtcDQDg61u8THJTeS9qblR3lSa8FRrM1mMh6fSIlMoSc9XVwyYYnNQDBgwAAmT55MWloaEyZMQKVSuaJcwkMp+x4X1wdxvqgjYDqDSMlF5DHbXTqTtTxGZuJqDxDHAQ5daLj2oLJGQ+bPo1nD4lYVDKx1PynjE0pAiGofZbJ6WloL7mEzEOTl5fHZZ5/xxhtvMH/+fO68805mzZrFTTfd5IryCQ+kbGqzb/Bii+ctbXcJrairyBpLXUhgCBBxlu45tZ1RbY6QsbMTjGz9Xa7mA9ZqtdpwTgaY3cdmIOjQoQMzZ85k5syZXLx4kUWLFnHzzTej1WpdUT7hoWqOHqUm4Cj+PXo0OGdtZzNLWUtbbWvBmLUAAYaxhjFXtrm2TG6mBATj1dMywOw+diWd++abb1i3bh2ff/45SUlJfPTRR84ul/Bgysyi+upqNDauNZ5ZZM6nWgvWJM7j0JfvUF2rbZC2Qq2+wr/ybKeyaI1jDMK17BojGDZsGDNnzuSvf/0rnTp1ckW5hAdTBo53PrjG5haXYNpCMOYzA8s2BHduz/mqq82611dnJAnHshkIfvzxR4KCLOdpF74tXHOCcr9eNre4tMbSwLIxn+g2AkIDAwgNDGDdvBtNjuu7TW60cpeeK7fYFK2X1UDwl7/8hd/97nf84Q9/sDhT6PXXX3dqwYTn66fJp9g/Enu2uLTGeGDZmNdnNm0qCyuXh6nVUNi10dueuXCJj2utT001J91IwhKrgUDJJ5SYmOiywgjvU19ZwbnCcptbXFpjrduoVWU2tcXSymU7Daw/yd3tYDnTbF5rvNK5MRIsfI/VQHDnnXcCkJCQwPXXX++yAgnvEZSaSujKLzkXGMSZgjaGNNWtLjOps1mZVbTfjv0I2q26gzho0K1kibLSuTHWgoUEh9bN5hjBk08+SUlJCTNmzGDWrFnExVmcDS18ULe0mQzcuJE++1+ldMgUDrQf3ezxAuF8ljbYMWcpWNjbkjDXnODR2CY7stjMeWwGgm3btnH27Fk++ugj0tPTqaioIC0tjT/+8Y+uKJ/wcMpU0p67VtN1WDDnCqOa3U1krtWluPYCloKFPS0Jc40FD2sBorFNdmSxmXPZtY6gV69ePP7444wfP56//OUvLF++XAKBAExzECndROeL2gEt6yJSUlz7+tRSu1jaSAcabKbTXPa0JMxZCx7mAcJ0rUQEU4cts/hestjMuWwGgiNHjrBu3To+/vhjevTowaxZs3jllVdcUTbhRYy7ifYNW4zmXMMVx02hTC2VVoEN1gaazRLgNfmZLQwg1oJHY60LWRPhPnZlH73nnnv48ssvCQsLc0WZhJdqyopje0irwA428hs1WRP2bG4O8wBhvFZC1kS4T6OBQKvVct1117Fo0aJmPTwrK4tFixah1WpZsGABS5ZYnhK4Z88eRo0axbp167j77uZPpRPuZbzi2BHMF5z5ygIzh2gsv1FjLHUxiVav0UDg5+fHhQsXqK2tpV27dk16sFar5dFHH+Wrr74iIiKCpKQkpkyZQmxsbIPrnn76adn+spWp0HXhw2sBIVxzgvjJ0RY3vLFFWXBmvAOaWq2m9Jssk2skQDiQtTGHpnLQGIVwPptdQ/369WPMmDFMmTLFJM/Qk08+2eh9OTk5REVFERkZCcCsWbPIzMxsEAj+8Y9/MH36dPbs2dOc8gsPFBXbiWOHLwFQ7teLcr9e1K5cw8CNGw3XBKWm2hUYzHdAMyeJ6xysBYvbTDi5i0k4ls1AEBYWRlhYGPX19VRWVtr94OLiYvr06WN4HRERwe7duxtcs2HDBrZu3SqBoBVJemIqSdf+fOi7YrJX55E3aDbntGcB/WpkNpTTZssaomI7kfTEVJvPNF6BbJy6WBLXOVhzu5TMSReTV7EZCJYuXdqsB+t0ugbHzHMWLV68mJdeegk/P79Gn5WRkUFGRgYARUVFZGdnN7k8VVVVzbrP23hiPXsnqrh0SkcNvQBo06EDbSorqdIFkb+7iPrnnuPK2LF2P8+kjv4BdA6LQK1We1y9W8oTP0t7DVOr6VxVSNXfxzQ4Vxo6jpIwfVewcR3V6isAFuusbGDjrT8PT/8sbQaC8ePHW0w6t3Xr1kbvi4iI4PTp04bXRUVFDWYd5ebmMmvWLADOnz/P5s2b8ff351e/+pXJdenp6aSnpwP63Ee2lt1bkm3Hcv3WwCPrmWz58PolX1DmH03l1lwG5hc0OG+t+8i8jqXfZHHuZCGl32S1qvECj/ws7dV5ARz4uGE6wlPb6XrpIIOu/gTov+C7dtVf1VN7iR0dxpOc3HCnNmU7S2/9eXj6Z2kzELz88suGP9fU1PDf//4Xf3/b69CSkpIoKCigsLCQ8PBw1q5dy5o1prNJCgsLDX9+4IEHSE1NbRAEROsVe8cQysy6jRT11dX0ycpnXJrt51gaUDY+11oCg1dpxrTW/nUnrm3Q03CF8cl2FWj9i712YZlarTbZm7m5BncfzNMjnnZAiUzZ/EYfPny4yesxY8Zw8803236wvz8rVqxg0qRJaLVa5s+fT1xcHCtXrgRg4cKFzSyyaC2Ulcf5OaWYp7Iuyy+juL6jXc+xNqBsKTCABAe3MgsQxon1ql6/hY5WNuipLk+gY3dXFNA32QwE5eXlhj/X19ezd+9ezp4928gdv0hJSSElxTS6WwsA//73v+16pmhd4saGW0xF8WEz1iKYp7S2NNPIODhIQPAs1jboAfT7LdTexqrJtrOseiKv7xoaPnw4KpUKnU6Hv78/AwYM4J133nFF2YRoEUt7HSjBQQKCEL+wGQiM+/GFcCV79kNuKvNuJJ/bCU0IC6wGgj179tCnTx969dJP+Xv//ff573//S79+/Xj22Wfp3l067ITzKPshH9500Cn7GygBYd2yJQ32TJYWgmc6XFLhtfmITLOsNl9sWBBL73T8njBtrJ146KGHDGklvv32W5YsWcLcuXPp0qWLYSqnEM4SPzmarup8ytRtOfRd03LhN0XMmGR69h9geH3uZKHFFczCvaYOCye2d5C7i9FqWW0RaLVaw2/969atIz09nenTpzN9+nSGDRvmqvIJH9UtbSZ9sp5HTbTTWgXQcBzBp/ZK9iLN2RPBkxhnWfVEVlsEWq0WjUafTPjrr79mwoQJhnPKcSGcSWkVaC5ccOn7Kl1FP23Jsn2xEK2A1RbBPffcw80330xwcDAdOnRg7LUUAMeOHaNLly4uK6DwXd3SZtJmi2NSWturscVpxtfIGIJoTawGgj/84Q/ccsstlJSUcNtttxnSTNTX1/OPf/zDZQUUQklp7d+jB/49e6JW13Nx7z7DeUfsj6ywN9upPeMIEjCEt2h0+uioUaMaHIuOjnZaYYQwp6S0rq+soLayAs2FC/hrNNRcW9RY7teLMwVqqnNz7cpiai9LaxDA8iI1S2RaqvAmdm1eL4S7KCmtL677iIqNG6HGNFHZKf9oDrQfTf7uEkLmzG1wv737HtjLWoAwJ4POwptIIBBeQdkGE6AwO5uh15br9wNKl3xBua4/O+tNpxfWV1YQuvJLhl+7XwhhmQQC4fVi7xhyLXFdiMnxMwVq1F2jKd2QbzLorIw1OHJsQQhvJoFAeD1riesOfVfM4U0H0Wh/aSkoYw3lF3XUHD1K54wsh3cfCeFtJBCIVstSgFDGGnbWB1HRpjvb68fZve9BUxmnrpAZRMKTSSAQPkUZa6j6rpj8nFLK8uF0FZyaM9ehLQNlPQLYnnIqQUK4mwQC4ZOU1sL6JV9QTgTbL42GDeW02bKGqNhOLZ6Kajy7qLEppzLNVHgCCQTCpykDzZpzKjQXLlCh68Kxw5dIcuB7NDbldN2yJRQdPshPW7IkGACcPQCr7rB8Lv5uy9tfihaTQCB8mvk4QnN2RmuJmDHJFB0+yFdvrWiw1zL+AS4ti9vF32393NkD+v9LIHAKCQRCuJHSCjAOAkp3UejNPtZCsLbhPVhvJQiHkEAghJmKNt3Z8Mo+l60zsJQKu+jwQdr07A0evM+taD2spqEWwheFa07QuaqIsvwyDm866JYyKDOOyguOuOX9he+RFoEQRuInR9Nv40a2XxpNuTbC4piBI2YVNSbh1skc2ZFNybECu3IW+cz008YGkl2hFQ9WSyAQwoiyzqDs75kcO1zR4Hy5Xy9y8uDYtQARrjlBP02+xWe1ZF1CzJhk1Gq1zet8Zt+ExgaSXeHUdv1/Bz5u1u3D1Goo7NrycvSKh9tfbPlzzEggEMICJeupuT1/z+TY4cuAPiiU+/Wi2D+ywXX11dUtWrGccOtkyv0DSLYxRuAz+yY0NpDsCrmrmh0EvIEEAiGawDhAHLq2Ohm6NriuLL+M4vqOTi9PS/dNsDdgeEWwcKYWBqL92dk2g7o7SSAQopmsJbsD169HMGfvvgn2BAxZ/dz6SSAQwofZEzDWLVtikkDPXj7fivAiEgiEcJL66mpOXds1zZtTXRsn0LOXzwxitxISCIRwAv8ePdAA1EDN0aOA9+6SZm83kzFrXU7SzeSZJBAI4QT+PXuivtqBfYMXUxNwlJCinWBhT2Xw7taCNdaCh+zl7JlkZbEQThA9IpTgiM4AVLYPoSxitMXrao4epWLjRlcWTYgGpEUghBMYzyja8Mo+zhe1Y9/gxQ2uqwk4ql+U5uLyCWFMAoEQThY9ItTquYo23eVfoXA7+SsohJPZWm9Q7teLb+c9b5KqIig1FUJDXFVE4eMkEAjhRlGxncjJgwPtRxtSVSjpKbhfAoFwDQkEQrhR0hNT6WiWqkJJT9HbrSUTvsSpgSArK4tFixah1WpZsGABS5aYTh1bvXo1L730EgCdO3fmX//6F0OHDnVmkYTwOJa2y6zQdUHz3s98uNp2qgpnp8UWrZ/Tpo9qtVoeffRRPv/8cw4fPsyHH37I4cOHTa4ZMGAA33zzDT/99BN/+tOfSE9Pd1ZxhPAaUbGdCFJdsutafVrsQPb8PdPJpRKtmdNaBDk5OURFRREZqe/3nDVrFpmZmcTGxhquGT36l7nVo0aNoqioyFnFEcJrKBlOs+3IWLnn75nk5AWSkxdo2CPBHsb7KLTGBW2iaZwWCIqLi+nTp4/hdUREBLt377Z6/TvvvMPtt9/urOII0SolPTEVjPZIsIfxPgr1lRWErvySgTYWtUmwaN2cFgh0Ol2DYyqVyuK127Zt45133mH79u0Wz2dkZJCRkQFAUVER2dnZTS5PVVVVs+7zNr5QT1+oIzShntd3off1Xex+bvtjOi6d0lFDL676B1PfoSM9z1vfdKVdQQHVe/aQn5/HlbFj7X4fS5Rd15R6yWfpGZwWCCIiIjh9+rThdVFREWFhYQ2u++mnn1iwYAGff/45PXr0sPis9PR0w/hBYmJiszZ4sKeZ3Rr4Qj19oY7gxHoaPVK/6jmAE0OeNRyLHhFqMnh9cd1HnF26lKDVa+iVX9Citz5VV4F/jx6Gesln6RmcFgiSkpIoKCigsLCQ8PBw1q5dy5o1pn2YP//8M9OmTeM///kP0dHRziqKEMIK81XPZwrUnClQX5vOqohC86tXCS3dY3V/ZnvVV1frs7IKj+K0QODv78+KFSuYNGkSWq2W+fPnExcXx8qVKwFYuHAhy5cv58KFCzzyyCOGe3Jzc51VJCGEGfOpq79sv2mqTN2WsvajuTAkpUXvV7V3BVrtBdYtW6Lf58A/oEXPE47h1HUEKSkppKSY/sVZuHCh4c9vv/02b7/9tjOLIIRoAmvpMKwFiKZStYvFr/awYeOazmERlH6TJZvVuJmsLBZC2NRYvqSm+PDBo9A2jLi07hzZkY1arTbZzUwCgntIIBBCuJyycU12djbdNTUc2ZEtu5e5kWxMI4Rwq4RbJ5O29EV69h/g7qL4LGkRCCE8xrmThR67nWVr7raSQCCE8AgxY5LdXQSrWnu3lQQCIYRHsLbhvSfw1FaKo0ggEEK4VEWb7mx4ZV+jW3gK15JAIIRwmXDNCfCHMwXtOFOgpmNPuLh3n+G8eXoLT9KS8Qu1Wk3pN1ktLkNIv0jGP+D4dP0SCIQQLtNPk0/owU8pHTKFYv9INEUaas6eBfQtBc25cx4ZCDx5/MIRJBAIIVwmKDUV0AeEfpp81Go1Xbt2BWB7/Tg09R3dWDrrWjp+4bNJ54QQwly3tJkm+xoUZmcz9NoX5M4mbKwjHEsCgRDCY9RXV3Nqzlx3F8Oi1rw5jwQCIYRH8O/Rg/I23dlZ73lTSOurq+mTlc+4NHeXxDkkEAghPELsHUOuZTgNcXdRGijLL+N0Fc1urXRTqzn1zrstLkf7mMH0+v3vW/wccxIIhBAewVEZTp1h/ZIvKCei2a0VTbAGf/+Wf912v6ijV4uf0pAEAiGEsKGlrRW1Wk3na7OjWqJjn84tfoYlEgiEEMKGlrZW9NNHb3BgiRxL0lALIYSPk0AghBA+TgKBEEL4OAkEQgjh4yQQCCGEj5NAIIQQPk4CgRBC+DgJBEII4eNUOp1O5+5CNEVwcDD9+/dv8n3nzp2jZ8+eji+Qh/GFevpCHcE36ukLdQTPqOfJkyc5f/68xXNeFwiaKzExkdzcXHcXw+l8oZ6+UEfwjXr6Qh3B8+spXUNCCOHjJBAIIYSP85lAkJ6e7u4iuIQv1NMX6gi+UU9fqCN4fj19ZoxACCGEZT7TIhBCCGGZTwSCrKwsBg0aRFRUFC+++KK7i9Mi/fv3Jz4+nmHDhpGYmAhAeXk5EydOZODAgUycOJGLFy8arn/hhReIiopi0KBBfPHFF+4qtk3z588nJCSEIUOGGI41p1579+4lPj6eqKgoHn/8cTypwWupjs8++yzh4eEMGzaMYcOGsXnzZsM5b6zj6dOnGT9+PDExMcTFxfHaa68Bre+ztFZPr/08da2cRqPRRUZG6o4fP667evWqLiEhQXfo0CF3F6vZ+vXrpzt37pzJsd/+9re6F154QafT6XQvvPCC7ne/+51Op9PpDh06pEtISNDV1NToTpw4oYuMjNRpNBqXl9ke33zzjW7v3r26uLg4w7Hm1CspKUm3c+dOXX19vW7y5Mm6zZs3u74yVliq49KlS3V//etfG1zrrXU8c+aMbu/evTqdTqerqKjQDRw4UHfo0KFW91laq6e3fp6tvkWQk5NDVFQUkZGRtGvXjlmzZpGZmenuYjlUZmYm999/PwD3338///d//2c4PmvWLNq3b8+AAQOIiooiJyfHjSW1bty4cXTv3t3kWFPrVVJSQkVFBTfeeCMqlYq5c+ca7vEElupojbfWsXfv3txwg34nrsDAQGJiYiguLm51n6W1elrj6fVs9YGguLiYPn36GF5HREQ0+oF5OpVKxW233cbw4cPJyMgAoLS0lN69ewP6v6BlZWWA99e9qfUqLi4mIiKiwXFPt2LFChISEpg/f76hy6Q11PHkyZP88MMPjBw5slV/lsb1BO/8PFt9INBZ6G9TqVRuKIlj7Nixg3379vH555/zxhtv8O2331q9trXVXWGtXt5Y34cffpjjx4+zf/9+evfuzVNPPQV4fx2rqqqYPn06r776KkFBQVava2319NbPs9UHgoiICE6fPm14XVRURFhYmBtL1DJK2UNCQrjrrrvIyckhNDSUkpISAEpKSggJCQG8v+5NrVdERARFRUUNjnuy0NBQ/Pz8aNOmDQ8++KCh686b61hXV8f06dO59957mTZtGtA6P0tr9fTGz7PVB4KkpCQKCgooLCyktraWtWvXMmXKFHcXq1kuX75MZWWl4c9ffvklQ4YMYcqUKbz33nsAvPfee0ydOhWAKVOmsHbtWq5evUphYSEFBQWMGDHCbeVvqqbWq3fv3gQGBvL999+j0+l4//33Dfd4KuXLEWDDhg2GGUXeWkedTsevf/1rYmJiePLJJw3HW9tnaa2eXvt5unx42g02bdqkGzhwoC4yMlL3/PPPu7s4zXb8+HFdQkKCLiEhQRcbG2uoy/nz53UTJkzQRUVF6SZMmKC7cOGC4Z7nn39eFxkZqYuOjvaoWRfmZs2apevVq5fO399fFx4ernv77bebVa89e/bo4uLidJGRkbpHH31UV19f747qWGSpjvfdd59uyJAhuvj4eN2dd96pO3PmjOF6b6zjd999pwN08fHxuqFDh+qGDh2q27RpU6v7LK3V01s/T1lZLIQQPq7Vdw0JIYRonAQCIYTwcRIIhBDCx0kgEEIIHyeBQAghfJwEAtGqXLhwwZD5sVevXoZMkJ07d+aRRx5x+PutXLmS999/v9n3P/DAA3z88ccOLJEQTefv7gII4Ug9evRg//79gD4lcOfOnfmf//kfp73fwoULnfZsIVxFWgTCJ2RnZ5OamgroA8T999/PbbfdRv/+/fnkk0/43e9+R3x8PJMnT6aurg7Q54m/+eabGT58OJMmTTJZNap49tlnefnllwFITk7m6aefZsSIEURHR/Pdd981uF6n0/HYY48RGxvLHXfcYUi+BrB8+XKSkpIYMmQI6enp6HQ6jh8/bshyCVBQUMDw4cMBWLJkCbGxsSQkJDg12InWTwKB8EnHjx9n06ZNZGZmct999zF+/HgOHDhAhw4d2LRpE3V1dfzmN7/h448/Zu/evcyfP58//OEPNp+r0WjIycnh1VdfZdmyZQ3Ob9iwgby8PA4cOMBbb73Fzp07Decee+wx9uzZw8GDB7ly5QobN27kuuuuo0uXLoZWzqpVq3jggQcoLy9nw4YNHDp0iJ9++ok//vGPDvvZCN8jgUD4pNtvv522bdsSHx+PVqtl8uTJAMTHx3Py5Eny8vI4ePAgEydOZNiwYTz//PMmycGsUZKPDR8+nJMnTzY4/+2333LPPffg5+dHWFgYEyZMMJzbtm0bI0eOJD4+nq1bt3Lo0CEAFixYwKpVq9Bqtaxbt47Zs2cTFBREQEAACxYs4JNPPqFjx44O+KkIXyVjBMIntW/fHoA2bdrQtm1bQ+rfNm3aoNFo0Ol0xMXFsWvXrmY918/PD41GY/EaS2mGa2pqeOSRR8jNzaVPnz48++yz1NTUADB9+nSWLVvGhAkTGD58OD169AD0my59/fXXrF27lhUrVrB169YmlVUIhbQIhLBg0KBBnDt3zhAI6urqDL+ht8S4ceNYu3YtWq2WkpIStm3bBmD40g8ODqaqqspkJlFAQACTJk3i4YcfZt68eYA+D/6lS5dISUnh1VdfNXQdCdEc0iIQwoJ27drx8ccf8/jjj3Pp0iU0Gg2LFy8mLi6uRc+966672Lp1K/Hx8URHR3PzzTcD0LVrVx588EHi4+Pp378/SUlJJvfde++9fPLJJ9x2220AVFZWMnXqVGpqatDpdPz9739vUbmEb5Pso0J4gZdffplLly7x3HPPubsoohWSFoEQHu6uu+7i+PHjMgYgnEZaBEII4eNksFgIIXycBAIhhPBxEgiEEMLHSSAQQggfJ4FACCF8nAQCIYTwcf8ftk4M2r5XGAcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv = rsf.predict_survival_function(X_test_sel, return_array=True)\n",
    "\n",
    "for i, s in enumerate(surv):\n",
    "    plt.step(rsf.unique_times_, s, where=\"post\", label=str(i))\n",
    "plt.ylabel(\"Survival probability\")\n",
    "plt.xlabel(\"Time in days\")\n",
    "plt.legend()\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can also plot the predicted cumulative hazard function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAELCAYAAADURYGZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6w0lEQVR4nO3dfVyUdb74/9dwLwISIAiOKYQ3gCBxo9Wad5vlsqWZ3didlRVZWnaqU31392R29lftbm121k6F27rVydQ8edxFdLdSNlcNxJtVpBQNTFARhBEQERiu3x/jTMMww3Az9/N+Ph7nsTFzzTXvT9P5vK/PvUpRFAUhhBBey8fZAQghhHAuSQRCCOHlJBEIIYSXk0QghBBeThKBEEJ4OUkEQgjh5fycHUBfRUVFMWrUqD5/7sKFCwwePNj2AbkYbyinN5QRvKOc3lBGcI1yVlZWUldXZ/Y9t0sEo0aNoqSkpM+fKywsZNq0abYPyMV4Qzm9oYzgHeX0hjKCa5QzKyvL4nvSNSSEEF5OEoEQQng5SQRCCOHl3G6MwJz29naqqqpobW21eM2QIUP49ttvHRhVd0FBQajVavz9/Z0ahxBCGPOIRFBVVUVoaCijRo1CpVKZvaapqYnQ0FAHR/YjRVE4d+4cVVVVxMfHOy0OIYQw5RFdQ62trURGRlpMAq5ApVIRGRnZY6tFCCGcwSMSAeDSSUDPHWIUQngfj+gachVbt25l6dKlaLVaHnnkEV588UVnhySEsIHDO6o5WlzT789rNJ007N034DiiRoRw/Z1jBnwfUx7TInA2rVbL4sWL2bJlC2VlZXz66aeUlZU5OywhhA0cLa6hrqrZ2WHYjbQIbKS4uJjExEQSEhIAmD9/Pps2bSI5OdnJkQnhHgb61G1PdVXNRKlDmPtsRr8+r1tZ3L/POoLdEsHJkydZsGABZ86cwcfHh9zcXJYuXdrlGkVRWLp0KQUFBQQHB/PnP/+ZjAzX/ZfVk+rqakaMGGH4W61WU1RU5MSIhHAv+qfuKHWIs0PpJjzwIpGlX3Pi/hX9+vwVGg0nPvjTgOMITBrHsF/8YsD3MWW3RODn58ebb75JRkYGTU1NZGZmMnPmzC5PyFu2bKG8vJzy8nKKiop4/PHHB1x5Lv/rYcpONXZ7XavV4uvr2697JseFseyWlB6vMXf0swwOC9E3A3nqtqcT9y+g9bvvYNw4Z4diF3ZLBLGxscTGxgIQGhpKUlIS1dXVXRLBpk2bWLBgASqVimuuuQaNRsPp06cNn3MnarWakydPGv6uqqoiLi7OiREJ4V46amvpOHeu30/d9tT63XcEjRvHyI8/6tfnKwoLmeDCm+s5ZIygsrKS/fv3M2nSpC6vm+tOqa6uHlAisPTkbu8FZdnZ2ZSXl1NRUcHw4cNZu3Yta9assdv3CeFpOs6do7OlxSWnsASNG0fYzTc7Owy7sXsiaG5uZt68eaxYsYKwsLAu7/W2OyUvL4+8vDxA96RdWFjY5f0hQ4bQ1NTUYxxardbqNQP129/+lpkzZ6LVarn//vu58soru31na2trt/htqbm52a73dwXeUEbwjnIal7GjowMCAqh4YKFzg+pJP38Pl/8tFTtqa2tTbrzxRuXNN980+35ubq6yZs0aw99jxoxRTp061eM9MzMzu71WVlZmNZbGxkar1zhCb2IdiO3bt9v1/q7AG8qoKN5RTuMyrnnkE2XNI584Lxg7coXf0lzdqWe3RpiiKDz88MMkJSXxzDPPmL1m9uzZfPTRRyiKwjfffMOQIUPccnxACNF/DevWc+L+BbpuIeEUdusa2rlzJx9//DGpqamkp6cD8Oqrr/LDDz8AsGjRInJycigoKCAxMZHg4GBWr15tr3CEEC6gYd16GvPzDX9fodFwprwcAJ/rZ+EXGems0Lya3RLB5MmTzY4BGFOpVLzzzjv2CkEI4WIa8/Np/e47asbPptovgY6oDvyG+eEXGUnzpUFEDXW9NQTeQFYWCyEc5oTfGKrTZ1HvOwyA4GEQEh4OQBQwZmKM84LzYpIIhBAOU+2XQKNPBHGjwxkzMYZabblLb73gLVxwxq4QwpOFddYz99kMUq4f7uxQxGWSCGxk4cKFREdHM378eGeHIoQQfSJdQzby4IMPsmTJEhYsWODsUIQQNnbwy618u7Ow35/XaDTU/GPrgOOIHpnA9AdzB3wfU9IisJEpU6YQERHh7DCEEHbw7c5CaisrnB2G3UiLQAghemHoqHjuWvZ6vz6rO49gmm0DsiHPSwRbXoQzh7q9PEjbAb79LO6wVPhZ//4DEEIIVyddQ0II4eU8r0Vg4cn9op23oRZCDMxAB2TtqbaygqGj4p0dht1Ii8BG7r77bq699lqOHDmCWq3mgw8+cHZIQrgVVx6QHToqnqSfTHN2GHbjeS0CJ/n000+dHYIQbm8gA7Ki/6RFIIQQXk5aBEIIp9KPDXh6P7wrk0QghHA4feWv0WhoPlUFgDp5vEf3w7sySQRCCIe50HaMi+2VfLHqLAAhcWpDAki7YZaTo/NekgiEEA5zsb2Sdm2DofKv9wty6RW33kIGi4UQDuXvewV3LXtdWgAuRBKBjZw8eZLp06eTlJRESkoKb7/9trNDEkKIXpGuIRvx8/PjzTffJCMjg6amJjIzM5k5cybJycnODk0IIXokicBGYmNjiY2NBSA0NJSkpCSqq6slEQiPcHhHNUeLawZ8nw788KPDBhEJW5JEYAeVlZXs37+fSZMmOTsUIfrEUoV/qlwDQNzo8AHd348OApWLA7qHsD2PSwS/Kf4N39V/1+11rVaLr69vv+45LmIcL0x8oVfXNjc3M2/ePFasWEFYWFi/vk8IZzlaXENdVTNR6pAur+sPmx/oOcMfzq8f0OeFfXhcInCm9vZ25s2bx7333sttt93m7HCE6JcodQhzn81wdhjCgTwuEVh6cm+y8zbUiqLw8MMPk5SUxDPPPGO37xHClky7gsy1BoTnk+mjNrJz504+/vhjtm3bRnp6Ounp6RQUFDg7LCF6pO8K0otShzBmYowTIxLO4HEtAmeZPHkyiqI4Owwh+sxaV1DDuvU05ufb5Ls6W1rwCQ62yb2E7UgiEEIYmKv0W/bsASA4O3vA9/cJDsYvMnLA9xG2JYlACC9jPC5QW1FP6KWznLh/BWC+0g/Ozibs5pu54q47B/zdQctfHPA9hO1JIhDCA/W0AMx4TUDopbNEV+2CKN17tqz0hfuQRCCEG+vPArDo8HZiavYw8sBRWr/7jqBx4xj58Ud2jlS4MkkEQrixo8U1hu4dYxHA8I7vGXngaLfP6Lt/yM4maNw4wm6+2QGRClcmiUAIF3R4RzUVX3XSsHdfj9fVVtQzuO4Y6Qfe7vVgrnT/CFOSCGyktbWVKVOmcOnSJTo6Orj99ttZvny5s8MSbujwjmoKPzkCQEh7LR3nzlm8dnBTIzE1JQxbvlwqdtFvkghsJDAwkG3bthESEkJ7ezuTJ0/mZz/7Gddcc42zQxNupmxzKeBPYsOXXFm4Eehh6qY/hC2Sp3sxMJIIbESlUhESolua397eTnt7OyqVyslRCXdgOuBb36AQ3vw9w1pKpRtHOIQkAhvSarVkZmZy7NgxFi9eLNtQC6uMu4H0M3zCOusZHnCGhmefYYKc5yscwOMSwZlXX+XSt923oe7Qaqnv5zbUgUnjGPaLX1i9ztfXlwMHDqDRaJg7dy6lpaWMHz++X98pvIO+JTDt3rGGLZ71i7sqmOyssISX8bhE4ArCw8OZNm0aW7dulUQgrIoObyck7/9xIk/3t35uvxCO4nGJwNKTu723oa6trcXf35/w8HAuXrzIl19+yQsv9O4wG+G9OmpraausoOXAHsOAsMztF45mt0SwcOFC8vPziY6OprS0tNv7hYWFzJkzh/j4eABuu+02XnrpJXuFY3enT5/mgQceQKvV0tnZyZ133snN8v/Mwgr91FCz0z8LCx0fkPBKdksEDz74IEuWLGHBggUWr7n++uvJt9H2ts6WlpbG/v37nR2GcEE97fvT6BNBWCgyK0g4ld0OppkyZQoRERH2ur0QbsP08BdjYZ31DO/43sERCdGVU8cIdu/ezYQJE4iLi+ONN94gJSXFmeEIYTeWDn/RzxASwpksJoLPP/+8xw8O9HD2jIwMTpw4QUhICAUFBdx6662Ul5ebvTYvL4+8PN2UiqqqKgpN+k6HDBlCU1NTj9+n1WqtXuMIra2t3eK3pebmZrve3xW4QxnrjymcP6E7sa5VA0HhmI35Co0GgAoz77lDOftKc7m8+nJ5YhnNcfVyWkwEf/3rXwE4e/Ysu3btYsaMGQBs376dadOmDTgRhIWFGf45JyeHJ554grq6OqKiorpdm5ubS25uLgBZWVlMM1lk8+2331qdEWTvWUO9FRQUxNVXX223+xcWFnb79+Np3KGMG/fuo6P58kHw4TBmYoxhnYCxEx/8CcDswjF3KKclB7/cyrc7C7u93q5pYOioeEO53LmMfeHq5bSYCFavXg3AzTffTFlZGbGxsYBudszixYsH/MVnzpwhJiYGlUpFcXExnZ2dRMoRdsKDWDsL2F1YqtR7UlWmmymoTu66jmboqHiSfjLNRpEJW7E6RlBZWWlIAgAxMTEcPdp9j3NTd999N4WFhdTV1aFWq1m+fDnt7e0ALFq0iA0bNvDuu+/i5+fHoEGDWLt2rezNIzxST4e/u9riMXOVvqVKvSfq5PEk/WQaaTfMsmV4wk6sJoJp06Zx0003cffdd6NSqVi7di3Tp0+3euNPP/20x/eXLFnCkiVLeh+pm9BqtWRlZTF8+HCPmRorBqYxP99ihe/sxWOmFb+5Sl8qdc9nNRGsXLmSjRs38vXXXwO6/vq5c+faPTB39fbbb5OUlERjY6OzQxEuxNnHQVrq3jGt+KXS9049JoLOzk7S0tIoLS2Vyr8Xqqqq2Lx5M7/85S/5/e9/7+xwhBN11OoOlDlx/wq7d//0pg/fUveOVPwCrCQCHx8fJkyYwA8//MCVV17pqJjc1tNPP81vf/tbl5imKpyr49w5OltawMe+3T8Hv9zKF6tWAj334UuFL3pitWvo9OnTpKSkMHHiRAYPHmx4/S9/+YtdA+uvHeuPUney+ypOrVaLbz+3oY4aEcL1d47p8Rr9vkqZmZkuPV9YOI5PcDAjVw2sO8ja077+SX/mo0ukkhf9ZjURLFu2zBFxuL2dO3fyl7/8hYKCAlpbW2lsbOS+++7jf/7nf5wdmnBDB7/cypHN/0fzqSrA8tO+POkLW7CaCKZOneqIOGzG0pO7vReUvfbaa7z22muAbvHIG2+8IUlA9Itpd49U9MLerG46980335CdnU1ISAgBAQH4+vp2WRUshLAtfVfQlVNnctey1yUJCLuz2iJYsmQJa9eu5Y477qCkpISPPvrI4p5AQmfatGkuvZxc2IbV7aU76/t9b3XyeIYmT+j354Xoi17tPpqYmGgYbH3ooYe47rrr7B2XEC5Pv710lDqk23t92V7adEC4trKCoaPibRWmEFZZTQTBwcG0tbWRnp7O888/T2xsLBcuXHBEbEK4HONWgD4JDGR7aXPTP/X78fS/PSFE31hNBB9//DGdnZ2sXLmSt956i5MnT/K///u/johNCKex1O1zqlwDQNzocMIDLxJZ+rXZSr+3i8j0LQFz0z9lGrJwFKuJ4Ouvv+bWW28lLCzMMJU0Pz+fxMREuwcnhLNY6vaJDm8npmYPIw8cpWXPHt2Llw+dN9aXRWTq5PEyICycymoiePLJJ3nzzTf59NNPSUpKAuCll16Sg9mFR9K3BGor6gm9dJaMA1u7vG9c+QdnZxN2881y3rBwe1YTQXx8PB988AG33347L7/8MnfccQeKojgiNiEcrmxzKfUNCiHnKxlaUwLqru/3p/LvaXWwDAwLV2A1EahUKjIyMvjHP/7B3XffTVFREVqt1hGxuZ1Ro0YRGhqKr68vfn5+lJSUODsk0Ucd584R0tLCZP9dhC2yzdP+tzsLLVb4clCLcAVWE4H+UJqoqCj+9re/8cILL1BaWmr3wNzV9u3bzR63KdyHrfcI0ieBu5a9boPohFOUrIZDG/r98XSNBirCBx7HsFT4me3/O7K6snjz5s0/Xuzjw+9+9zs6OzttHogQnkTfCgB56vcIhzbAmUPOjsJurLYIamtr+c1vfkNZWRmtra2G17dt22bXwNyRSqXixhtvRKVS8dhjj5Gbm+vskIQTSSvAhgb4RD5gZw7pnsYf2mz9WjMOuOvh9Xr33nsvd911F5s3b+a9997jww8/ZOjQoY6IrV+2/zmPsye6r+jUdmjx9evfNtTRIxOY/qD1Sn3nzp3ExcVx9uxZZs6cybhx45gyZUq/vlMIYUT/RD4s1TnfPywVUm93znc7gNVEcO7cOR5++GHefvttpk6davg/0V1cXBwA0dHRzJ07l+LiYkkEQtjKAJ7IRc+sJgJ/f39AN2i8efNm4uLiqKqqsntg/WXpyd3e21BfuHCBzs5OQkNDuXDhAn//+9956aWX7PZ9wjZMVxAPdLM4IdyR1UTwq1/9ivPnz/Pmm2/y5JNP0tjYyFtvveWI2NxKTU2N4Vznjo4O7rnnHmbNktWirs50BXFfNosTwlNYTQT6FcRDhgxh+/btdg/IXSUkJPCvf/3L2WGIfohShzBNfYzG/Hy7HzQvetDTgLAzxwe8QK9mDa1atYrKyko6OjoMr//pT3+ya2BCOEJHbS0d585xZtV/Aj+uHBZO0NOAsIcP1jqb1UQwZ84crr/+em644YZ+H/4uhKvRjw3UNyiEtLTIvkGuQgaEncJqImhpaeE3v/mNI2IRwmH0YwNhnfUMDzjDyNV9X0ksewgJT9GrMYKCggJycnIcEU+/KYqCSqVydhg9ks36HM90VpC+K0g/OyjjwIo+jQkYV/5VZbqtVvQHyhiT1cTCnVhMBKGhoahUKhRF4dVXXyUwMBB/f39DhdvY2OjIOHsUFBTEuXPniIyMdNlkoCgK586dIygoyNmheI3DO6op/OQIoDtIBnSbynW2tBAWDMM7vrd4boClp33jyl+dPJ6kn0yTswSE27OYCJqamhwZx4Co1Wqqqqqora21eE1ra6vTK+GgoCDUarX1C0WfWDtNLPXSLkYeOAr8eHKYpU3l9AnA0tO+VP7CE/Xq8HpX5+/vT3x8z/2xhYWFXH311Q6KSDiS/gwB04Vg4U2NxNSUMFLdZnjN2slh+s3ipMIX3sQjEoHwXod3VHNW40948/dc5/N11zf96deZArJZnPA2kgiES9N3++gHeUG3cvvTT9YAUO87DIAR/Zz5I4ToxXkEAP/85z9ZvXo1oFtgVlFRYdeghNAr21zK2aNnaausoLOp+wSFCO0ZUi/tInXWGCdEJ4RnsNoiWL58OSUlJRw5coSHHnqI9vZ27rvvPnbu3OmI+ISX63J05OUFX4U23NvddHaQzP8X3shqIti4cSP79+8nIyMD0G217E4zioT7s8XRkZaYnics8/+FN7KaCAICAlCpVIb5+RcuXLB7UEI4kgwOC29ndYzgzjvv5LHHHkOj0bBq1SpuuOEGHn30UUfEJoQQwgGstgiee+45vvjiC8LCwjhy5AivvPIKM2fOdERswsvteWsT9b7DiNCesel9jccFZExAiF4kgrfeeos77rhDKn/hcMfKLoBvKInJg212z4NfbuWLVSsB3SphGRMQoheJoLGxkZtuuomIiAjmz5/P7bffTkxMjCNiE16sYd16OpsaiQiF7H+7p1ef6Wk3UD391hEzH10iq4aFuMxqIli2bBnLli3j4MGDrFu3jqlTp6JWq/nyyy97/NzChQvJz88nOjqa0tLSbu8risLSpUspKCggODiYP//5z4aZSUIc2noUTfh1RIe3W73W2v5AxmTrCPe0pugHNh2odnYY/abRXOTdI7sHfJ/kuDCW3ZJig4i66vXK4ujoaIYNG0ZkZCRnz561ev2DDz7IkiVLWLBggdn3t2zZQnl5OeXl5RQVFfH4449TVFTU+8iFR6v2SwAg+eeWK3U92R/I8206UE3Z6UaSY8OcHYpHspoI3n33XdatW0dtbS233347q1atIjk52eqNp0yZQmVlpcX3N23axIIFC1CpVFxzzTVoNBpOnz5NbGxsnwogPFeE9gwp18/o1bUyBdRNmJxLnK7RQEW47g8r5xInx4ax7rFr7RufnegWQbpu7FYTwYkTJ1ixYgXp6ek2/eLq6mpGjBhh+FutVlNdXS2JQAhP1s9ziRt8v+a8bzEPbXXPFoFGo+HDrR8O+D7jIsbxwsQXbBBRVxYTQWNjI2FhYTz//PMA1Nd33eI3IiJiQF9s7rQuS4fK5OXlkZeXB0BVVRWFhYV9/r7m5uZ+fc7duHM5648pnD+h++/ikuoKBrfXmi2LcRlry/5FVVkpIXFqty23Je78W1qSrtFA0AgOxP87oCtjSEjIjxc0A2bKfE7ZRTun0Gjc8zwPrVaLRqMZ8H2qWqoobCkc8H1MWUwE99xzD/n5+WRmZhpOKtNTqVR8//33A/pitVrNyZMnDX9XVVURFxdn9trc3Fxyc3MByMrK6tc+M7bcn8aVuVs5jQ+VOX35IJm40eG0nvmO4fxAeEd0t5lAGo2G8PBw4MdZQNf+/FbS3KjcveFuv2UXJl1ABq0nYViqoVzGZexpQLhNgWC/K9k4/3/tFLB9ufpvaTER5OfnA9htp9HZs2ezcuVK5s+fT1FREUOGDJFuIS9jfJRkhPYMEeiOjxx54KjhJLFvdjb3uOhLBohdlKUuoB66f3oaEA4O9CNqcKA9IhX0Yozgpz/9KV999ZXV10zdfffdFBYWUldXh1qtZvny5bS366YCLlq0iJycHAoKCkhMTCQ4ONiwzbXwHmWbSwF/xh5Zw2ijU8TA6CSx7/Z1Gwh29acrt2bpSb6v9Engoc19+pilAWF3HRtwFxYTQWtrKy0tLdTV1dHQ0GDoGmpsbOTUqVNWb/zpp5/2+L5KpeKdd97pY7jCHVk6U7i+QSH8/FEyF91o+RSx5fvsHJ3ooqfB3L7o4clfuB6LieD9999nxYoVnDp1iszMTEMiCAsLY/HixQ4LULg34+6fuNHhAIbTxkJaWhgxqI4r7lrkxAhFN/14khfuzWIiWLp0KUuXLuUPf/gDTz75pCNjEm7OuAVw6vIAcOqlXYw8cBSAlj17AAjOziZsluWD5A9+uZWqstIeVwoLC/rbxWOL1oBwO1bHCJ588klKS0spKyujtbXV8LqlFcNClG0upb5BIayznvCmRmJqShhpNAYQnJ1tOG3MEuPN4WRTuB5YqvBP/FP3vyMn9+1+0qXjlXp1VGVhYSFlZWXk5OSwZcsWJk+eLIlAWKTv9rnO52vwh7BFPVf65uinjMrmcFZY6tMfOVlXoWc95Jy4hFuxmgg2bNjAv/71L66++mpWr15NTU0NjzzyiCNiE27MFsdLqpPHSxLoDRfs0+/tJnHGm7HJXkLOYzURDBo0CB8fH/z8/GhsbCQ6OnrAi8mEZ2pYt57G/Hw6O6fgExzs7HCEE/Vnk7jk2DCuHHWQh7bmdXvvSP0RxkaMtWWIwojVRJCVlYVGo+HRRx8lMzOTkJAQJk6c6IjYhJtpzM+n9bvv8EmfhV9kpLPDEf1gq+2e9UnA2iZxxpuxfXb0M17ZrRsXyorJ6nLd2Iix5CTkDDguYZ7VRPDf//3fgG4R2KxZs2hsbCQtLc3ugQn30rBuPeVVAdSmP01ziJqooSHWPyRcjq22e06ODWNO+vA+fabg+wIAXrr2Je4Yc8eAvl/0jcVEsG+f5YU8+/btk0NkBPBjd1DLnj3UpC/lQmA0Q9UhjJnYu1PsLJ0q5lVnCVuY+dNli2ZLejnds7dP+r19kreXrJgsSQJOYDERPPvssxY/pFKp2LZtm10CEu5F3x1Ue+29aALHEBcfztxne/+QoD9UxrTS96qzhAewmrdm8Gg2aTL46v2eT78qqtDtHjwpvuddg/vzJC/cn8VEsH37dkfGIdyQcXdQve8wgF63BKDrgjGvP1TGzMyfA73YU+mp93frnuKtjM1Pio9gTvpw7pl05QADFZ7I6hjBRx+ZnwIo6whEY34+NTHXcSEwmrj4cMZMjCHl+u5Pk5a6f/RbSHvNk78FNU2t1DVf4hWTp/renHPr7K4c4RmsJoI9l7cDAN1GdF999RUZGRmSCAQn/MagGdxzd5DxCmHTrSJkC2mgZDUx9SX8gPXjX82RrhxhC1YTwR/+8Icuf58/f57777/fbgEJ96E/YN60O8i4BaB/6vfqFcI97ftzeSuI/UNu6PZU7+rn3NrKZ0c/o+D7Alkr4ERWE4Gp4OBgysvL7RGLcEPmDpg3HgD22qd+48q/p31/Rk4mT5PBV8E55DouOqfRV/p6Go2GYyeOAboZQ7JWwDmsJoJbbrnFcJZwZ2cnZWVl3Hln3/aNEZ6lNyuITQ+T8QqWKn+TfX9Mp3KWnbc+2OtuTCt8vZKaEqDrgjF9ApBpo85jNRE899xzP17s58fIkSNRq93zAGlhG5ZWEOu7hLxmDYBpl4+Fyn9N0Q9s2lsNe3UDv6ZTOT2xn99SV49ppS+nzbkGq4lg6tSpgO5kso6ODgDq6+uJiOh5PrLwbEHjxhE0bhyNZ0tYt/xF4MfxAH13kEfqqcvH3JP/+7u7VfzuOJXT0hO+JfoksHqWHEHrDqwmgry8PP7jP/7DsPmcoiioVCrZeM5LGdYOJM6guaqZtqaDNJ2t9ezxgL50+Vh48nelir+vlTqY79LpiewN5F6sJoLf/e53HD58mKioKEfEI1yIfizAmPFWEoEB39HYVOlZC8LMzfDpofLXW1P0A7/YeAhw3Sd/fQLoa6Wuv1b68T2X1URw1VVXESxbCnsN48rf+EhJveDsbAKGxxM8NIK25i8BD1sQZm67hx4qf/2gr/7p/9W5qS5T8RvT7ez5CiCVuujOaiJ47bXXuO6665g0aRKBgYGG1//rv/7LroEJ59APBAeNG2f2SMnDO6o5+8kR4obq/vaYw2P0LQF9ErBy0ItpC8BZT/+97ebRtwJkZ09hjtVE8NhjjzFjxgxSU1Px8fFxREzCyYLGjWPkx123FtEfSK8/jH7MxBhKPWnfQeMk0Isze/UtAXu3AKxV9L3t5pFWgOiJ1UTg5+fH73//e0fEIpxI3yWkbw2YOlpcQ11VM3GjwxkUfITSbX93/2mixuMBvWwJGJsUH2GTJGCpsjddbGWOVPDCFqwmgunTp5OXl8ctt9zSpWtIpo96FuMkEHbzzcCPrQCAuqpmotQhzH02g3XL1xuSgFuPDxi3AnpoCZjby9+W5+v2tL2CVPTCEawmgjVr1gC6sQI9mT7qmYLGjaM59zX2FdfAm/sM3UBxo8OJMjlsxu1XDpes1s0GGjm5SyvAXKVvbi9/Wy8CMzfnXhZbCUexmggqKiocEYdwooZ162nZs4fg7GxDF1CUOoS40d23ljY+Q8At6buD9FNCL7cC9AnAXKXvatNAhbA1OY9A0JifT3XsT2gYfgcaoy4gc/S7irpNl1BP20AYTQnVn9XriErfdExAdt0UzibnEQgA3UrhS4O6dQEZM24NuMWU0ZLVkP+07p8tbANhzFEHvJiOCcgqXOFsch6B4ITfGOp9hxFn0hIwPVnM7U4U07cEbl5htuKHH7uEbDn4q2dpNpDswyNcjZxH4KWMVxCf7JwCgbq1AeYOldGPB7jlXkIjJ/eYBIwXhQ108Ne04rc0x19aAMLVyHkEXsp4uqhPcDDR4e2kXD+cdcv/4FqHypjZ+yddo4GKcOufNd0qwoStF4WZdvnI1E/hLuQ8Ai+mX0G87819hu2k9UnAZaaGmtv7p7d6sUq4v4vCzHX7SJePcFcWE8GxY8eoqakxnEegt2PHDi5dusRVV11l9+CE7em7hI7VhXFWfR373txHncl20i4zBmBhrv+BPsyv158JYM5AxgXMLQKTLh/hriwmgqeffppXX3212+uDBg3i6aef5q9//atdAxP2cWjrUU52TkEzUnfwfBwQpQ6hoSqAweEu1BKAH7uEerH3jznmtoY21p9FYaYHrcvTv/AEFhNBZWUlaWlp3V7PysqisrLSnjEJO2lYt56TF6NoHqLusljs4Jdb+aL4KIPDXXCRWA+DvdbYagzAuBvIeABYnv6Fp7CYCFpbWy1+6OLFi3YJRtiXbpbQdURcoeoyTdQlF4ld7haqicjiKZOuHY3mIu8eMd/dY0y/QGygA8HGLQAZABaeyGIiyM7OZtWqVTz66KNdXv/ggw/IzMy0e2DC9k74jUEzeAxxQ8O7vee0RWLmTgQDwwrgTdrr+t2Xb8v9gKQbSHgyi4lgxYoVzJ07l08++cRQ8ZeUlNDW1sbGjRsdFqAYGP3g8Am/MRwKvA7A4sphh+hpywdjIydTFDKDV/eOY1J81xW/us3YbLsCuKd9/2ULCOHpLCaCmJgYdu3axfbt2ykt1S0s+vnPf86MGTMcFpwYOMPgcKBucHji2CbDuIC+S8ih5wqYTgftYcuH37+/G6i36S6fpnpzjq/MBhKerlfnEUyfPr1fN9+6dStLly5Fq9XyyCOP8OKLL3Z5v7CwkDlz5hAfr6uEbrvtNl566aV+fZcwr9ovgeaAiG47iX67s9CQABw+ZbQPB8DY6vAXc+QcXyF0+rzFRG9ptVoWL17MF198gVqtJjs7m9mzZ5OcnNzluuuvv578y1sdCNtqWLeezqZGwkJh7rO6/n99S8BhC8dMu4L6uzjMDvRdQXKOr/B2dksExcXFJCYmkpCg65KYP38+mzZt6pYIhP3oZwn5RUYaXjNOAg5pBZh2BZms9jV3EIyePTaCM5UVkyVJQHg9uyWC6upqRowYYfhbrVZTVFTU7brdu3czYcIE4uLieOONN0hJSbFXSF7HdJaQ8TbSNm8JWJr9Y+YsYOPVvuYOgtGzxawfGQQWwjq7JQJFUbq9pt+8Ti8jI4MTJ04QEhJCQUEBt956q9mdTfPy8sjLywOgqqqKwsLCPsfT3Nzcr8+5G305B+3YwcmLQyEQzjft4P1/y6P5VBUAPkNjbfLvIvbU34ip+RqA8PO6CQWaISaL0oJGUBOYxunL31d4sp0/H24DYOwVPoy9wodr4/yYNuKS+S+5+D2FhV2PRe3Lb7nmzBqq26oZHtA9ocT4xJDYnuiy/114w3+z3lBGcP1y2i0RqNVqTp48afi7qqqKuLi4LteEhf3Y7M/JyeGJJ56grq6OqKioLtfl5uaSm5sL6FY29+ccV284/7Vh3Xp++J//ITw8nJY9eyB9KcE++zn37XbAhttImx73OHIyhOtm/4Sbmf0TDuifu9+9PBNoIKt9e/tbfnb0M46dOEZWTJZbrgHwhv9mvaGM4PrltFsiyM7Opry8nIqKCoYPH87atWtZs2ZNl2vOnDlDTEwMKpWK4uJiOjs7iTTqzxZ905ifj19VFYSHU3vtvWgCx0Db5wDMfHSJbRaMmZ76ZWHqp57pGEBfV/ua69rRaDR8uPVD66FenhIqUz+F6JndEoGfnx8rV67kpptuQqvVsnDhQlJSUnjvvfcAWLRoERs2bODdd9/Fz8+PQYMGsXbt2m7dR6L3TviN4UT6DYSMH8+pcg0Ag0IDiBrRz1XD5vr99a2AHk79MmZ6+ldf+/3N7fLZWzIlVIjesVsiAF13T05O16exRYsWGf55yZIlLFmyxJ4heJVqvwQuqK4gBIgbHc6g4CMcLuzjZnLGlb+5Vb8WWgGWZv/ok0B/zgL+7OhnlNSUdOvacfVmthDuxq6JQDiOfs3A4EEdzH32ZwCsW74esLKZXE9bPvSi60fP0rm/fWkBWDrqUbp2hLAvSQQeQr9moDM0FOg6VdRst5C5AV/9//ay8jduBfTlyd/SlE7TbR6ka0cIx5BE4MaMD6A/VhdG3bBWaPwb65YXGQ6e79YaMJcAelnxQ9fK33gNQF+e/C31+0vFL4RzSCJwUw3r1rP3vb9TE3MdPqFh1I8chrZpPSpqgXDzU0X7MOPHUp+/ceU/KT6COenD+zQN1FK/vxDCeSQRuKlDW49yZOw9gG5gOA5oqAqmnaHdVw2btgJMZvyYq/QtrfjtT+VvTN8lJP3+QrgOSQRuqtpPt4fTtHvHdjluMiRO3f1i/X4/ZloBls71HWiF3xPZ30cI1yKJwA3pZwgFBh6jdNvfKd2GYUwgYnTSjxfqWwJm9vvRs9W5vkII9yWJwA0d2nqUukEhdLQUc75Mt3WEfkyg3i/I6EKjJHB5x8+BrvTtL/1MIdnoTQjXI4nADVX7JaBt0W32Zrp1RLeNrYxaAua6gWx5rq855k4Ak/EBIVyLJAI3YXz2cH3gdfgrbUT38sB5fStAPwDsyG4gfStApoYK4bokEbgJ47OHOy4dpEN7Foju1Wf1q37tOQBsaZGYvitIpooK4bokEbiJo52dNCnFBHQcoKOlEujF1hEn/klNRBZFp+qZFB/Rr/1+rLF2+Lsc/C6E65NE4AYa1q2npfUYCk1EjRgLWDhXoGQ16fv/CBXhhjUDK2omANhlHEAOfxfCM0gicAOHth6lwy+YAL9BPR8xeWgDIc0VEH41jJxMniaDT2uu6deYQE9HPOrpWwFy+LsQ7k0SgRs42tmJ0lFFcPQYq9c2h8QTfnmW0Ffv72ZSPP0aE+jNVE9pBQjhGSQRuDh9txBA1s9vtHzh5TEBLp8ZvKboB4oq6s0eCm+JcStABnmF8B6SCFyYfmO5jrhgAv0GWZ4qarSZ3J9brmH3+7sNU0WtjQ0YV/7GA74yyCuE95BE4MIObT3K4VHjUVq+tNwtZJQElpPLLr8ZhNP7vYKMu4Ckq0cI7ySJwAXpF48d6biSjs5vAAvdQiZJoCz2Nv7f2EtMm9a3aaLSBSSEd5NE4GL2vLWJo0X1wHU0++wDum8jAZhNAnPSh8PF73u8v+lsINn7RwghicCF7HlrE8VHQiE8lOjwdgLqShkUOqbHJJA35CnKgnMMi8UKC7snAkvjACALvoQQkghcyrGyC+AbyojwYi4O0qBoaxkcHt/9wsuHzecNeYo/nJ9McrDle5ou+pJxACGEKUkELibw4j8pbygGMGwtbaqmqZWKziRerbmGSfGWdw81TgKy6EsIYYkkAhfRsG49zY2HaES3TbTZcQGAktXE1JdQQZLVFcP67iBJAkKInkgicDL9DKGWPXu4OCEb6J4E1hT9QPOuVfzk4nZS2nSJ4uAVM8ntxYphORZSCGGNJAIn0Ff+AC179gBQe+29dLRV4RMQwv93fAgc3w3AT1sKSGv4gmt8vgXgcEAqOwdNJ+S6R7vcUz8grNFo+HDrh4DMCBJC9I4kAgcyfvoHCM7OJjg7m7MZ89h/sAKlo4qLgSNIPv05twfoEkFK2yHwgZqILGKuu4+UrIdIMXNv/cKwGJ8Yw2syI0gI0RuSCBykYd16zixbBsCZjFROXRGKX/RQLmja0JRsQemoAmBINCwhD9qAkZOByZB6OzFZD5m9r+lZwA8EPcC0adMcUyghhEeQROAgjfn5HIy/mtroUC61VkNdC0GXBnHpYgcAV8SNIXKYijkX/6j7wM0rwELlb8w4CeQk5MApOxZCCOGRJBE4yB6GUBVWAa2NBIWOIiQyjbBo3aKuMRNjaD6zlkmHdVM9i1JeYpKVJGDaEtBvEVF4qtCu5RBCeB5JBA7QsG49p9rOA5Ay7V5mPX634b2iz96EHRuZdHk2UFHKS0y641mL9zJ3NKSMAwghBkISgQPoZwgFBg3vkgQAQso3MqLtOIcDUmkePZcfJlzJe1sttwZME4BMDRVCDJQkAjswnh4K8M15Fe2DzhHkPwrQrQvYdKAagOfatJwMuIqUX/yz23YQ5kgCEELYmiQCGzOeHRScrVsgdjbEF7TgM+pq7np/Nwk/fMa/+e4iNMiPA6FnKAi/goCtD8kZwEIIp5BEMECmT//fHf+OiqRxdIRFo/L3B6Ct+TQ+wVfyw8XD/Fv1Kq7x1y0OI3Yyb6haqfCFscjTvhDCOSQRDFBjfj6t331H0LhxAJwYeRUXlDZUHSq0im5qKAFDaQsIJiNiDe8PHsz7QaNh8FAIjeZIfYMcDCOEcCpJBP2kbwl8c17F2dHj8QnU7QXd2n4ale9QCqNv5WCglknxEfy0pYArOj/glahIALJiUg33kdW/Qghnk0TQD/pxgIPxV1MV1ghaCGIUAEEhsVSHjSVk2J/4SeQxAnx92K2cp2SQLglI/78QwtVIIuijbkmArmsD1hT9QNVfVhAYtp8jqgDGMgiChpAVPJScCQslCQghXI4kgh6YDgSDbrfQ6tifUB3hBx2NpEy7l68SGlm2eh4AQa21JIw4zZGAAMaGXsnqO7Y4I3QhhOg1H3vefOvWrYwdO5bExERef/31bu8risJTTz1FYmIiaWlp7Nu3z57h9Jl+INjY4ewbOBjnh6KtxS8intU+oyj+/hMU5Sgj24+ToDoNwNjQK8mZsNAZYQshRJ/YrUWg1WpZvHgxX3zxBWq1muzsbGbPnk1ycrLhmi1btlBeXk55eTlFRUU8/vjjFBUV2SukXvvLi69z4oejKFoVKqOBYIDWpgrd/w4OoHroXjqUr2kP7CClrY3VgWN0F6Xe3qsN44QQwhXYLREUFxeTmJhIQkICAPPnz2fTpk1dEsGmTZtYsGABKpWKa665Bo1Gw+nTp4mNjbVXWBbpK3+ANu1ZAAJUkSho0bY2G67z94+iLvoSX44/iLrDnyBfX8b6+JOTcCPc+JbD4xZCiIGyWyKorq5mxIgRhr/VanW3p31z11RXV9slERxa9TEH3rM8V1+rnAPAVxWJryoCZVA7hZMO6F7zUXW59oS/FrVfLOsXfmnzOIUQwtHslggURen2mkql6vM1AHl5eeTl5QFQVVVFYWFhfwLq8W1fVSTaQR3svOYAABrVEOpVowgLUBEe2DWmOCBjcFb/4rCz5uZml4zLlryhjOAd5fSGMoLrl9NuiUCtVnPy5EnD31VVVcTFxfX5GoDc3Fxyc3MByMrK6vcJXN5wcldhYaHHl9MbygjeUU5vKCO4fjntNmsoOzub8vJyKioqaGtrY+3atcyePbvLNbNnz+ajjz5CURS++eYbhgwZ4pTxASGE8GZ2axH4+fmxcuVKbrrpJrRaLQsXLiQlJYX33nsPgEWLFpGTk0NBQQGJiYkEBwezerXstyOEEI5m1wVlOTk55OR03Udn0aJFhn9WqVS888479gxBCCGEFXZdUCaEEML1SSIQQggvJ4lACCG8nCQCIYTwcpIIhBDCy6kUc8t7XVhUVBSjRo3q8+dqa2sZOnSo7QNyMd5QTm8oI3hHOb2hjOAa5aysrKSurs7se26XCPorKyuLkpISZ4dhd95QTm8oI3hHOb2hjOD65ZSuISGE8HKSCIQQwst5TSLQb1rn6byhnN5QRvCOcnpDGcH1y+k1YwRCCCHM85oWgRBCCPO8IhFs3bqVsWPHkpiYyOuvv+7scAZk1KhRpKamkp6eTlZWFgD19fXMnDmT0aNHM3PmTBoaGgzXv/baayQmJjJ27Fj+9re/OStsqxYuXEh0dDTjx483vNafcu3du5fU1FQSExN56qmnzB5+5Czmyvjyyy8zfPhw0tPTSU9Pp6CgwPCeO5bx5MmTTJ8+naSkJFJSUnj77bcBz/stLZXTbX9PxcN1dHQoCQkJyvHjx5VLly4paWlpyuHDh50dVr+NHDlSqa2t7fLav//7vyuvvfaaoiiK8tprrynPP/+8oiiKcvjwYSUtLU1pbW1Vvv/+eyUhIUHp6OhweMy98Y9//EPZu3evkpKSYnitP+XKzs5Wdu3apXR2diqzZs1SCgoKHF8YC8yVcdmyZcrvfve7bte6axlPnTql7N27V1EURWlsbFRGjx6tHD582ON+S0vldNff0+NbBMXFxSQmJpKQkEBAQADz589n06ZNzg7LpjZt2sQDDzwAwAMPPMD//d//GV6fP38+gYGBxMfHk5iYSHFxsRMjtWzKlClERER0ea2v5Tp9+jSNjY1ce+21qFQqFixYYPiMKzBXRkvctYyxsbFkZGQAEBoaSlJSEtXV1R73W1oqpyWuXk6PTwTV1dWMGDHC8Ldare7xB3N1KpWKG2+8kczMTMM5zjU1NYaT3WJjYzl79izg/mXva7mqq6tRq9XdXnd1K1euJC0tjYULFxq6TDyhjJWVlezfv59JkyZ59G9pXE5wz9/T4xOBYqa/TaVSmbnSPezcuZN9+/axZcsW3nnnHb7++muL13pa2fUslcsdy/v4449z/PhxDhw4QGxsLM8++yzg/mVsbm5m3rx5rFixgrCwMIvXeVo53fX39PhEoFarOXnypOHvqqoq4uLinBjRwOhjj46OZu7cuRQXFxMTE8Pp06cBOH36NNHR0YD7l72v5VKr1VRVVXV73ZXFxMTg6+uLj48Pjz76qKHrzp3L2N7ezrx587j33nu57bbbAM/8LS2V0x1/T49PBNnZ2ZSXl1NRUUFbWxtr165l9uzZzg6rXy5cuEBTU5Phn//+978zfvx4Zs+ezYcffgjAhx9+yJw5cwCYPXs2a9eu5dKlS1RUVFBeXs7EiROdFn9f9bVcsbGxhIaG8s0336AoCh999JHhM65KXzkCbNy40TCjyF3LqCgKDz/8MElJSTzzzDOG1z3tt7RUTrf9PR0+PO0EmzdvVkaPHq0kJCQov/71r50dTr8dP35cSUtLU9LS0pTk5GRDWerq6pQZM2YoiYmJyowZM5Rz584ZPvPrX/9aSUhIUMaMGeNSsy5MzZ8/Xxk2bJji5+enDB8+XPnjH//Yr3Lt2bNHSUlJURISEpTFixcrnZ2dziiOWebKeN999ynjx49XUlNTlVtuuUU5deqU4Xp3LOOOHTsUQElNTVUmTJigTJgwQdm8ebPH/ZaWyumuv6esLBZCCC/n8V1DQggheiaJQAghvJwkAiGE8HKSCIQQwstJIhBCCC8niUB4lHPnzhl2fhw2bJhhJ8iQkBCeeOIJm3/fe++9x0cffdTvzz/44INs2LDBhhEJ0Xd+zg5ACFuKjIzkwIEDgG5L4JCQEJ577jm7fd+iRYvsdm8hHEVaBMIrFBYWcvPNNwO6BPHAAw9w4403MmrUKD7//HOef/55UlNTmTVrFu3t7YBun/ipU6eSmZnJTTfd1GXVqN7LL7/MG2+8AcC0adN44YUXmDhxImPGjGHHjh3drlcUhSVLlpCcnMzPf/5zw+ZrAK+88grZ2dmMHz+e3NxcFEXh+PHjhl0uAcrLy8nMzATgxRdfJDk5mbS0NLsmO+H5JBEIr3T8+HE2b97Mpk2buO+++5g+fTqHDh1i0KBBbN68mfb2dp588kk2bNjA3r17WbhwIb/85S+t3rejo4Pi4mJWrFjB8uXLu72/ceNGjhw5wqFDh1i1ahW7du0yvLdkyRL27NlDaWkpFy9eJD8/n6uuuoohQ4YYWjmrV6/mwQcfpL6+no0bN3L48GEOHjzIr371K5v9uxHeRxKB8Eo/+9nP8Pf3JzU1Fa1Wy6xZswBITU2lsrKSI0eOUFpaysyZM0lPT+fXv/51l83BLNFvPpaZmUllZWW397/++mvuvvtufH19iYuLY8aMGYb3tm/fzqRJk0hNTWXbtm0cPnwYgEceeYTVq1ej1WpZt24d99xzD2FhYQQFBfHII4/w+eefExwcbIN/K8JbyRiB8EqBgYEA+Pj44O/vb9j618fHh46ODhRFISUlhd27d/frvr6+vnR0dJi9xtw2w62trTzxxBOUlJQwYsQIXn75ZVpbWwGYN28ey5cvZ8aMGWRmZhIZGQnoDl366quvWLt2LStXrmTbtm19ilUIPWkRCGHG2LFjqa2tNSSC9vZ2wxP6QEyZMoW1a9ei1Wo5ffo027dvBzBU+lFRUTQ3N3eZSRQUFMRNN93E448/zkMPPQTo9sE/f/48OTk5rFixwtB1JER/SItACDMCAgLYsGEDTz31FOfPn6ejo4Onn36alJSUAd137ty5bNu2jdTUVMaMGcPUqVMBCA8P59FHHyU1NZVRo0aRnZ3d5XP33nsvn3/+OTfeeCMATU1NzJkzh9bWVhRF4a233hpQXMK7ye6jQriBN954g/Pnz/Of//mfzg5FeCBpEQjh4ubOncvx48dlDEDYjbQIhBDCy8lgsRBCeDlJBEII4eUkEQghhJeTRCCEEF5OEoEQQng5SQRCCOHl/n+5I8cKVfA0YwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv = rsf.predict_cumulative_hazard_function(X_test_sel, return_array=True)\n",
    "\n",
    "for i, s in enumerate(surv):\n",
    "    plt.step(rsf.unique_times_, s, where=\"post\", label=str(i))\n",
    "plt.ylabel(\"Cumulative hazard\")\n",
    "plt.xlabel(\"Time in days\")\n",
    "plt.legend()\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutation-based Feature Importance\n",
    "\n",
    "The implementation is based on scikit-learn's Random Forest implementation and inherits many\n",
    "features, such as building trees in parallel. What's currently missing is feature importances\n",
    "via the `feature_importance_` attribute.\n",
    "This is due to the way scikit-learn's implementation computes importances. It relies on\n",
    "a measure of *impurity* for each child node, and defines importance as the amount of\n",
    "decrease in impurity due to a split. For traditional regression, impurity would be measured by the variance, but for survival analysis\n",
    "there is no per-node impurity measure due to censoring. Instead, one could use the\n",
    "magnitude of the log-rank test statistic as an importance measure, but scikit-learn's\n",
    "implementation doesn't seem to allow this.\n",
    "\n",
    "Fortunately, this is not a big concern though, as scikit-learn's definition\n",
    "of feature importance is non-standard and differs from what Leo Breiman\n",
    "[proposed in the original Random Forest paper](https://github.com/scikit-learn/scikit-learn/pull/8027#issuecomment-327595859).\n",
    "Instead, we can use permutation to estimate feature importance, which is\n",
    "preferred over scikit-learn's definition. This is implemented in the\n",
    "[permutation_importance](https://scikit-learn.org/stable/modules/generated/sklearn.inspection.permutation_importance.html#sklearn.inspection.permutation_importance)\n",
    "function of scikit-learn,\n",
    "which is fully compatible with scikit-survival."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.inspection import permutation_importance\n",
    "\n",
    "result = permutation_importance(rsf, X_test, y_test, n_repeats=15, random_state=random_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>importances_mean</th>\n",
       "      <th>importances_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>pnodes</th>\n",
       "      <td>0.077655</td>\n",
       "      <td>0.019356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>0.015749</td>\n",
       "      <td>0.008232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>progrec</th>\n",
       "      <td>0.009565</td>\n",
       "      <td>0.013541</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>horTh=yes</th>\n",
       "      <td>0.006289</td>\n",
       "      <td>0.004057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tgrade</th>\n",
       "      <td>0.003346</td>\n",
       "      <td>0.003764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>menostat=Post</th>\n",
       "      <td>0.000070</td>\n",
       "      <td>0.001023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tsize</th>\n",
       "      <td>-0.000070</td>\n",
       "      <td>0.009698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>estrec</th>\n",
       "      <td>-0.008220</td>\n",
       "      <td>0.009413</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               importances_mean  importances_std\n",
       "pnodes                 0.077655         0.019356\n",
       "age                    0.015749         0.008232\n",
       "progrec                0.009565         0.013541\n",
       "horTh=yes              0.006289         0.004057\n",
       "tgrade                 0.003346         0.003764\n",
       "menostat=Post          0.000070         0.001023\n",
       "tsize                 -0.000070         0.009698\n",
       "estrec                -0.008220         0.009413"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(\n",
    "    {\n",
    "        k: result[k]\n",
    "        for k in (\n",
    "            \"importances_mean\",\n",
    "            \"importances_std\",\n",
    "        )\n",
    "    },\n",
    "    index=X_test.columns,\n",
    ").sort_values(by=\"importances_mean\", ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result shows that the number of positive lymph nodes (`pnodes`) is by far the most important\n",
    "feature. If its relationship to survival time is removed (by random shuffling),\n",
    "the concordance index on the test data drops on average by 0.077655 points.\n",
    "Again, this agrees with the results from the original\n",
    "[Random Survival Forests paper](https://projecteuclid.org/euclid.aoas/1223908043)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
